本稿の目的

本稿では、博士論文の2つ目の論文となる群れ内の凝集性に関する分析を行う。

0. パッケージの読み込み

本稿はRの基本操作とtidyverseパッケージによるデータハンドリングができることを前提としている。tidyverseパッケージを用いたデータ処理については、以下の書籍などを参照。

使用するパッケージは以下のとおりである。

## データハンドリング
library(tidyverse)
library(readxl)
library(knitr)
library(easystats)
library(DT)
library(scales)
library(emmeans)
## モデリング  
library(brms)
library(rethinking)
library(rstan)
library(rstanarm)
library(cmdstanr)
library(ggeffects)
library(DHARMa)
library(DHARMa.helpers)
## ネットワーク分析  
library(sna)
library(ANTs)
library(asnipe)
library(aninet)
library(vegan)
library(igraph)
library(tidygraph)
library(ggraph)
## グラフや表関連
library(ggtext)
library(ggforce)
library(ggbeeswarm)
library(plotly)
library(bayesplot)
library(viridis)
library(ggnewscale)
library(GGally)
library(flextable)
library(ggrepel)
library(patchwork)
library(kableExtra)
library(ggsci)
library(lemon)
library(ggsignif)
library(ggh4x)
library(ggrepel)
library(ggpattern)
# ## フォント関連
library(extrafont)
require(systemfonts)
# ベイズモデルの設定
rstan_options(auto_write = TRUE) 
options(mc.cores = parallel::detectCores()) 

1 データの読み込み

以下、基本的なデータの読み込みを行う。なお、本論文で使用するのはた8つの期間(2018交尾期、2019非交尾期、2019交尾期、2020非交尾期、2020交尾期、2021非交尾期、2021交尾期、2022非交尾期)のデータである。

1.1 Daily data の読み込み

以下、日ベースのデータについてデータフレームを作成する。
なお、**には観察年の下2桁と交尾期(m)/非交尾期(nm)が入る。

  • base_**: 追跡時間や気温についての情報。
  • female_**: 各メスの確認の有無と発情の有無、年齢など(分派集団で確認された個体も含む)
  • male_**: 各オスの確認の有無
  • female_all: 全期間をまとめたメスの確認状況
  • male_all: 全期間をまとめたオスの確認状況

1.2 個体追跡データの読み込み

以下では、個体追跡の生データとして、以下のデータシートを読み込む。

  • focal_18m: 2018年交尾期の個体追跡データ
  • focal_19m: 2019年交尾期の個体追跡データ
  • focal_20m: 2020年交尾期の個体追跡データ
  • focal_21m: 2021年交尾期の個体追跡データ
  • focal_19nm: 2019年非交尾期の個体追跡データ
  • focal_21nm: 2021年非交尾期の個体追跡データ
  • focal_22nm: 2022年非交尾期の個体追跡データ
  • focal_combined: 全個体追跡データを結合したもの

1.3 群れの詳細な出欠データを読み込む

1日に2つ以上の集団を確認した日を含めて、群れ個体の詳細な出欠データを読み込む。
5歳以上のメスと群れオスの確認状況のデータが記録されている(2018年度は6歳以上のメス)。

以下のデータを読み込んで加工する。
なお、**には観察年の下2桁と交尾期(m)/非交尾期(nm)が入る。

  • group_**: 観察した各集団における個体の確認状況と観察時間
  • group_all: 全てのデータを結合

1.4 オスの攻撃データを読み込む

以下、全生起サンプリングによって収集した生データを読み込む。

  • agg_m**: 各交尾期のデータ
  • aggression_all: 全交尾期のデータ

2 メス同士のCSIの算出

2.1 CSIの定義

ここでは、TYとのCSIを算出する。

CSIの算出式は以下の通り。なお、\(i = 1,2,3,...,N\)はメスのIDを表す。\(P_i\)の算出の際には互いに毛づくろいしているポイントは除いた。対象としたのは2018年時点で6歳以上だった非発情の個体である。

  • \(G_{ij}\): メスiとメスjの毛づくろい時間割合 (iとjが毛づくろいしたポイント数/iとj総ポイント数)
  • \(P_{i}\): メスiとメスjの近接時間割合 (3m以内にいとjがいたポイント数/iまたはjが地上採食または休息したポイント数)
  • \(\overline{G}: \frac{1}{N} \sum^{N}_{i=1}G_{i}\)
  • \(\overline{P}: \frac{1}{N} \sum^{N}_{i=1}P_{i}\)

\[ CSI_{i} = \frac{\left(\frac{G_{i}}{\overline{G}} +\frac{P_{i}}{\overline{P}}\right)}{2} \;\;(i = 1,2,\dots,N) \]

使用するのは2018年交尾期、2019年交尾期、2019年非交尾期、2021年非交尾期、2022年非交尾期のデータである。

2.2 データの加工

交尾期のデータについて、追跡個体が発情しているかの列を追加し、発情していない日のデータのみを抽出する(focal_18m_bfocal_19m_b)。

focal_18m %>% 
  left_join(female_18m %>% 
              rename(subject = femaleID) %>% 
              select(date, subject, rs2) , by = c("date", "subject")) %>% 
  filter(rs2 == "0") -> focal_18m_b

focal_19m %>% 
  left_join(female_19m %>% 
              rename(subject = femaleID) %>% 
              select(date, subject, rs2) , by = c("date", "subject")) %>% 
  filter(rs2 == "0") -> focal_19m_b

2.3 CSIの算出(全期間)

算出するメスのIDは、全期間のデータがあった以下の通り。

femaleID <- sort(c("Kil","Mik","Koh","Aka","Ten","Ntr","Hen","Hot","Tot","Mei","Ako","Kol","Mal","Kit","Kun"))

2.3.1 毛づくろい時間割合の算出

focal_all <- bind_rows(focal_18m_b, focal_19m_b, focal_19nm, focal_21nm, focal_22nm)

## 毛づくろい相手を表す列を作成  
focal_all %>% 
  mutate(groom = ifelse(groomer == subject,groomee,groomer), 
         groom2 = ifelse(groomer2 == subject,groomee2,groomer2)) -> focal_all_b

各メスの追跡時間(総瞬間サンプリングポイント数)を算出する。

focal_all_b %>% 
  filter(subject %in% femaleID) %>% 
  group_by(subject) %>% 
  summarise(dur = n()) -> focal_duration

対象となるメスと毛づくろいしていたポイントを抽出する。

focal_all_b %>% 
  filter(activity == "G") %>% 
  filter(subject %in% femaleID) %>% 
  select(date, no_focal, time, subject, groom, groom2) %>% 
  pivot_longer(cols = c("groom","groom2"), 
               names_to = "groom", values_to = "ID") %>% 
  filter(ID %in% femaleID) -> groom_pairs

最後に、毛づくろい時間割合を算出してマトリックスにする。

groom_mat <- df.to.mat(groom_pairs, 
                       actor = "subject",
                       receiver = "ID",
                       sym = T,
                       tob = focal_duration$dur)

2.3.2 近接時間割合の算出

続いて、近接時間割合を算出する。なお、互いに毛づくろいをしていた時間は除く点は注意が必要である。

まず、地上採食・休息・毛づくろいのポイントのみを抽出する。

focal_all_b %>% 
  filter(activity %in% c("F","R","G") & T_G == "G") %>% 
  select(no_focal, subject, x0_1m:x1_3m, groom, groom2, study_period) %>% 
  replace_na(list(groom = "NA", groom2 = "NA")) -> focal_prox

続いて、上記のうち各メス同士が毛づくろいをしていなかったポイント数を算出する(= 分母)。

prox_denom <- matrix(nrow = length(femaleID), ncol = length(femaleID))

for(i in seq_along(femaleID)){
  for(j in seq_along(femaleID)){
    prox_denom[i,j] <- focal_prox %>% 
                          filter(subject == femaleID[i]) %>% 
                          filter(groom != femaleID[j] & groom2 != femaleID[j]) %>% 
                          nrow()
  }
}

## 転置行列と足す
prox_denom2 <- prox_denom + t(prox_denom)

## 対角成分は0にする  
diag(prox_denom2) <- 0
rownames(prox_denom2) <- femaleID
colnames(prox_denom2) <- femaleID

続いて、各メスが毛づくろいをしていなかったポイントのうち近接していたポイント数を算出する。

prox_numer <- matrix(nrow = length(femaleID), ncol = length(femaleID))

for(i in seq_along(femaleID)){
  for(j in seq_along(femaleID)){
    prox_numer[i,j] <- focal_prox %>% 
                          filter(subject == femaleID[i]) %>% 
                          filter(groom != femaleID[j] & groom2 != femaleID[j]) %>% 
                          filter(str_detect(x0_1m, femaleID[j])|str_detect(x1_3m, femaleID[j])) %>% 
                          nrow()
  }
}

## 転置行列と足す
prox_numer2 <- prox_numer + t(prox_numer)

## 対角成分は0にする  
rownames(prox_numer2) <- femaleID
colnames(prox_numer2) <- femaleID

最後に、これらを基に近接時間割合を算出する。

prox_mat <- prox_numer2/prox_denom2
diag(prox_mat) <- 0

2.3.3 CSIの算出

最後に、CSIを算出する。ここでは、aninetパッケージのdyadic_csi()関数を用いる。

library(aninet)
CSI_female <- dyadic_csi(list(groom_mat, prox_mat))
colnames(CSI_female) <- femaleID
rownames(CSI_female) <- femaleID

CSIを図示すると以下のようになる。

CSI_female %>% 
  data.frame() %>% 
  rownames_to_column(var = "ID1") %>% 
  pivot_longer(2:16, 
               names_to = "ID2",
               values_to = "CSI") %>% 
  ggplot(aes(x = ID1, y = ID2))+
  geom_tile(aes(fill = CSI))+
  scale_fill_gradient2(high = muted("blue"), low = muted("red"), mid = "white",
                       midpoint = 1)+
  theme(aspect.ratio = 1)+
  labs(x = "", y = "")

2.4 CSIの算出(2018交尾期~2019交尾期)

算出するメスのIDは、この期間のデータがあった以下の通り。

femaleID2 <- sort(c("Kil","Mik","Koh","Aka","Tam","Ten","Ntr","Hen","Hot","Tot","Mei","Ako","Kol","Mal","Kit","Kun"))

2.4.1 毛づくろい時間割合の算出

focal_fh <- bind_rows(focal_18m_b, focal_19m_b, focal_19nm)

## 毛づくろい相手を表す列を作成  
focal_fh %>% 
  mutate(groom = ifelse(groomer == subject,groomee,groomer), 
         groom2 = ifelse(groomer2 == subject,groomee2,groomer2)) -> focal_fh_b

各メスの追跡時間(総瞬間サンプリングポイント数)を算出する。

focal_fh_b %>% 
  filter(subject %in% femaleID2) %>% 
  group_by(subject) %>% 
  summarise(dur = n()) -> focal_duration_fh

対象となるメスと毛づくろいしていたポイントを抽出する。

focal_fh_b %>% 
  filter(activity == "G") %>% 
  filter(subject %in% femaleID2) %>% 
  select(date, no_focal, time, subject, groom, groom2) %>% 
  pivot_longer(cols = c("groom","groom2"), 
               names_to = "groom", values_to = "ID") %>% 
  filter(ID %in% femaleID2) -> groom_pairs_fh

最後に、毛づくろい時間割合を算出してマトリックスにする。

groom_mat_fh <- df.to.mat(groom_pairs_fh, 
                          actor = "subject",
                          receiver = "ID",
                          sym = T,
                          tob = focal_duration_fh$dur)

2.4.2 近接時間割合の算出

続いて、近接時間割合を算出する。なお、互いに毛づくろいをしていた時間は除く点は注意が必要である。

まず、地上採食・休息・毛づくろいのポイントのみを抽出する。

focal_fh_b %>% 
  filter(activity %in% c("F","R","G") & T_G == "G") %>% 
  select(no_focal, subject, x0_1m:x1_3m, groom, groom2, study_period) %>% 
  replace_na(list(groom = "NA", groom2 = "NA")) -> focal_prox_fh

続いて、上記のうち各メス同士が毛づくろいをしていなかったポイント数を算出する(= 分母)。

prox_denom_fh <- matrix(nrow = length(femaleID2), ncol = length(femaleID2))

for(i in seq_along(femaleID2)){
  for(j in seq_along(femaleID2)){
    prox_denom_fh[i,j] <- focal_prox_fh %>% 
                            filter(subject == femaleID2[i]) %>% 
                            filter(groom != femaleID2[j] & groom2 != femaleID2[j]) %>% 
                            nrow()
  }
}

## 転置行列と足す
prox_denom_fh2 <- prox_denom_fh + t(prox_denom_fh)

## 対角成分は0にする  
diag(prox_denom_fh2) <- 0
rownames(prox_denom_fh2) <- femaleID2
colnames(prox_denom_fh2) <- femaleID2

続いて、各メスが毛づくろいをしていなかったポイントのうち近接していたポイント数を算出する。

prox_numer_fh <- matrix(nrow = length(femaleID2), ncol = length(femaleID2))

for(i in seq_along(femaleID2)){
  for(j in seq_along(femaleID2)){
    prox_numer_fh[i,j] <- focal_prox_fh %>% 
                            filter(subject == femaleID2[i]) %>% 
                            filter(groom != femaleID2[j] & groom2 != femaleID2[j]) %>% 
                            filter(str_detect(x0_1m, femaleID2[j])|str_detect(x1_3m, femaleID2[j])) %>% 
                            nrow()
  }
}

## 転置行列と足す
prox_numer_fh2 <- prox_numer_fh + t(prox_numer_fh)

## 対角成分は0にする  
rownames(prox_numer_fh2) <- femaleID2
colnames(prox_numer_fh2) <- femaleID2

最後に、これらを基に近接時間割合を算出する。

prox_mat_fh <- prox_numer_fh2/prox_denom_fh2
diag(prox_mat_fh) <- 0

2.4.3 CSIの算出

最後に、CSIを算出する。ここでは、aninetパッケージのdyadic_csi()関数を用いる。

library(aninet)
CSI_fh_female <- dyadic_csi(list(groom_mat_fh, prox_mat_fh))
colnames(CSI_fh_female) <- femaleID2
rownames(CSI_fh_female) <- femaleID2

CSIを図示すると以下のようになる。

CSI_fh_female %>% 
  data.frame() %>% 
  rownames_to_column(var = "ID1") %>% 
  pivot_longer(2:16, 
               names_to = "ID2",
               values_to = "CSI") %>% 
  ggplot(aes(x = ID1, y = ID2))+
  geom_tile(aes(fill = CSI))+
  scale_fill_gradient2(high = muted("blue"), low = muted("red"), mid = "white",
                       midpoint = 1)+
  theme(aspect.ratio = 1)+
  labs(x = "", y = "")

2.5 CSIの算出(Tam生存前まで)

算出するメスのIDは、この期間のデータがあった以下の通り。

femaleID2 <- sort(c("Kil","Mik","Koh","Aka","Tam","Ten","Ntr","Hen","Hot","Tot","Mei","Ako","Kol","Mal","Kit","Kun"))

2.5.1 毛づくろい時間割合の算出

focal_withTam <- bind_rows(focal_18m_b, focal_19m_b, focal_19nm, focal_20m, focal_21nm) %>% 
  mutate(date = as_date(date)) %>% 
  filter(date <= "2021-03-14")

## 毛づくろい相手を表す列を作成  
focal_withTam %>% 
  mutate(groom = ifelse(groomer == subject,groomee,groomer), 
         groom2 = ifelse(groomer2 == subject,groomee2,groomer2)) -> focal_withTam_b

各メスの追跡時間(総瞬間サンプリングポイント数)を算出する。

focal_withTam_b %>% 
  filter(subject %in% femaleID2) %>% 
  group_by(subject) %>% 
  summarise(dur = n()) -> focal_duration_withTam

対象となるメスと毛づくろいしていたポイントを抽出する。

focal_withTam_b %>% 
  filter(activity == "G") %>% 
  filter(subject %in% femaleID2) %>% 
  select(date, no_focal, time, subject, groom, groom2) %>% 
  pivot_longer(cols = c("groom","groom2"), 
               names_to = "groom", values_to = "ID") %>% 
  filter(ID %in% femaleID2) -> groom_pairs_withTam

最後に、毛づくろい時間割合を算出してマトリックスにする。

groom_mat_withTam <- df.to.mat(groom_pairs_withTam, 
                          actor = "subject",
                          receiver = "ID",
                          sym = T,
                          tob = focal_duration_withTam$dur)

2.5.2 近接時間割合の算出

続いて、近接時間割合を算出する。なお、互いに毛づくろいをしていた時間は除く点は注意が必要である。

まず、地上採食・休息・毛づくろいのポイントのみを抽出する。

focal_withTam_b %>% 
  filter(activity %in% c("F","R","G") & T_G == "G") %>% 
  select(no_focal, subject, x0_1m:x1_3m, groom, groom2, study_period) %>% 
  replace_na(list(groom = "NA", groom2 = "NA")) -> focal_prox_withTam

続いて、上記のうち各メス同士が毛づくろいをしていなかったポイント数を算出する(= 分母)。

prox_denom_withTam <- matrix(nrow = length(femaleID2), ncol = length(femaleID2))

for(i in seq_along(femaleID2)){
  for(j in seq_along(femaleID2)){
    prox_denom_withTam[i,j] <- focal_prox_withTam %>% 
                            filter(subject == femaleID2[i]) %>% 
                            filter(groom != femaleID2[j] & groom2 != femaleID2[j]) %>% 
                            nrow()
  }
}

## 転置行列と足す
prox_denom_withTam2 <- prox_denom_withTam + t(prox_denom_withTam)

## 対角成分は0にする  
diag(prox_denom_withTam2) <- 0
rownames(prox_denom_withTam2) <- femaleID2
colnames(prox_denom_withTam2) <- femaleID2

続いて、各メスが毛づくろいをしていなかったポイントのうち近接していたポイント数を算出する。

prox_numer_withTam <- matrix(nrow = length(femaleID2), ncol = length(femaleID2))

for(i in seq_along(femaleID2)){
  for(j in seq_along(femaleID2)){
    prox_numer_withTam[i,j] <- focal_prox_withTam %>% 
                            filter(subject == femaleID2[i]) %>% 
                            filter(groom != femaleID2[j] & groom2 != femaleID2[j]) %>% 
                            filter(str_detect(x0_1m, femaleID2[j])|str_detect(x1_3m, femaleID2[j])) %>% 
                            nrow()
  }
}

## 転置行列と足す
prox_numer_withTam2 <- prox_numer_withTam + t(prox_numer_withTam)

## 対角成分は0にする  
rownames(prox_numer_withTam2) <- femaleID2
colnames(prox_numer_withTam2) <- femaleID2

最後に、これらを基に近接時間割合を算出する。

prox_mat_withTam <- prox_numer_withTam2/prox_denom_withTam2
diag(prox_mat_withTam) <- 0

2.5.3 CSIの算出

最後に、CSIを算出する。ここでは、aninetパッケージのdyadic_csi()関数を用いる。

library(aninet)
CSI_withTam_female <- dyadic_csi(list(groom_mat_withTam, prox_mat_withTam))
colnames(CSI_withTam_female) <- femaleID2
rownames(CSI_withTam_female) <- femaleID2

CSIを図示すると以下のようになる。

CSI_withTam_female %>% 
  data.frame() %>% 
  rownames_to_column(var = "ID1") %>% 
  pivot_longer(2:16, 
               names_to = "ID2",
               values_to = "CSI") %>% 
  ggplot(aes(x = ID1, y = ID2))+
  geom_tile(aes(fill = CSI))+
  scale_fill_gradient2(high = muted("blue"), low = muted("red"), mid = "white",
                       midpoint = 1)+
  theme(aspect.ratio = 1)+
  labs(x = "", y = "")

2.6 CSIの算出(調査期間ごと)

2.6.1 毛づくろい時間割合の算出

各メスの追跡時間(総瞬間サンプリングポイント数)を算出する。

focal_all_b %>% 
  filter(subject %in% femaleID) %>% 
  group_by(subject, study_period) %>% 
  summarise(dur = n()) -> focal_duration_sp

対象となるメスと毛づくろいしていたポイントを抽出する。

focal_all_b %>% 
  filter(activity == "G") %>% 
  filter(subject %in% femaleID) %>% 
  select(date, no_focal, time, subject, groom, groom2, study_period) %>% 
  pivot_longer(cols = c("groom","groom2"), 
               names_to = "groom", values_to = "ID") %>% 
  filter(ID %in% femaleID) -> groom_pairs

最後に、毛づくろい時間割合を算出してマトリックスにする。

groom_mat_sp <- list()
sp <- unique(groom_pairs$study_period)

for(i in c(1,2,3,5)){
groom_mat_sp[[i]] <- df.to.mat(groom_pairs %>% filter(study_period == sp[i]), 
                             actor = "subject",
                             receiver = "ID",
                             sym = T,
                             tob = focal_duration_sp %>% filter(study_period == sp[i]) %>%.$dur)
}

## 2021非交尾期はホタルがfemaleIDに含まれるメスと毛づくろいをしなかったので、別に求める。  
groom_mat_sp[[4]] <- df.to.mat(groom_pairs %>% filter(study_period == sp[4]), 
                               actor = "subject",
                               receiver = "ID",
                               sym = T,
                               tob = focal_duration_sp %>% 
                                 filter(study_period == sp[4]) %>% 
                                 filter(subject != "Hot") %>%
                                 .$dur)

## ホタルの行列を追加  
groom_mat_sp[[4]] %>% 
  data.frame() %>% 
  bind_rows(data.frame(row.names = "Hot")) %>% 
  mutate(Hot = 0L)-> groom_mat_sp4

## 行名と列名を並べ替え  
groom_mat_sp4[15,1:15] <- 0  
groom_mat_sp[[4]] <- groom_mat_sp4 %>% 
                       as.matrix()
groom_mat_sp[[4]] <- groom_mat_sp[[4]][sort(rownames(groom_mat_sp[[4]])),sort(colnames(groom_mat_sp[[4]]))] 

2.6.2 近接時間割合の算出

続いて、近接時間割合を算出する。なお、互いに毛づくろいをしていた時間は除く点は注意が必要である。

続いて、上記のうち各メス同士が毛づくろいをしていなかったポイント数を算出する(= 分母)。

mat <- matrix(nrow = length(femaleID), ncol = length(femaleID))
prox_denom_sp <- list(mat,mat,mat,mat,mat)
prox_denom_sp2 <- list(mat,mat,mat,mat,mat)

for(i in seq_along(sp)){
  for(j in seq_along(femaleID)){
    for(k in seq_along(femaleID))
    prox_denom_sp[[i]][j,k] <- focal_prox %>% 
                                  filter(study_period == sp[i]) %>% 
                                  filter(subject == femaleID[j]) %>% 
                                  filter(groom != femaleID[k] & groom2 != femaleID[k]) %>% 
                                  nrow()
    
    diag(prox_denom_sp[[i]]) <- 0
    rownames(prox_denom_sp[[i]]) <- femaleID
    colnames(prox_denom_sp[[i]]) <- femaleID
    prox_denom_sp2[[i]] <- prox_denom_sp[[i]] + t(prox_denom_sp[[i]])
  }
}

続いて、各メスが毛づくろいをしていなかったポイントのうち近接していたポイント数を算出する。

prox_numer_sp <- list(mat,mat,mat,mat,mat)
prox_numer_sp2 <- list(mat,mat,mat,mat,mat)

for(i in seq_along(sp)){
  for(j in seq_along(femaleID)){
    for(k in seq_along(femaleID))
    prox_numer_sp[[i]][j,k] <- focal_prox %>% 
                                  filter(study_period == sp[i]) %>% 
                                  filter(subject == femaleID[j]) %>% 
                                  filter(groom != femaleID[k] & groom2 != femaleID[k]) %>% 
                                  filter(str_detect(x0_1m, femaleID[k])|str_detect(x1_3m, femaleID[k])) %>% 
                                  nrow()
    
    diag(prox_numer_sp[[i]]) <- 0
    rownames(prox_numer_sp[[i]]) <- femaleID
    colnames(prox_numer_sp[[i]]) <- femaleID
    prox_numer_sp2[[i]] <- prox_numer_sp[[i]] + t(prox_numer_sp[[i]])
  }
}

最後に、これらを基に近接時間割合を算出する。

prox_mat_sp <- list(mat,mat,mat,mat,mat)

for(i in seq_along(sp)){
  prox_mat_sp[[i]] <- prox_numer_sp2[[i]]/prox_denom_sp2[[i]]
  diag(prox_mat_sp[[i]]) <- 0
}

2.6.3 CSIの算出

最後に、CSIを算出する。ここでは、aninetパッケージのdyadic_csi()関数を用いる。

library(aninet)
CSI_sp_female <- list(mat,mat,mat,mat,mat)

for(i in seq_along(sp)){
  CSI_sp_female[[i]] <- dyadic_csi(list(groom_mat_sp[[i]], prox_mat_sp[[i]]))
  colnames(CSI_sp_female[[i]]) <- femaleID
  rownames(CSI_sp_female[[i]]) <- femaleID
}

CSIを図示すると以下のようになる。

CSI_sp_df <- data.frame()

for(i in seq_along(sp)){
CSI_df <- CSI_sp_female[[i]] %>% 
  data.frame() %>% 
  rownames_to_column(var = "ID1") %>% 
  pivot_longer(2:16,
               names_to = "ID2",
               values_to = "CSI") %>% 
  mutate(study_period = sp[i])

CSI_sp_df <- bind_rows(CSI_sp_df, CSI_df)
}

CSI_sp_df %>% 
  ggplot(aes(x = ID1, y = ID2))+
  geom_tile(aes(fill = CSI))+
  scale_fill_gradient2(high = muted("blue"), low = muted("red"), mid = "white",
                       midpoint = 1)+
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = -45,
                                   size = 7),
        axis.text.y = element_text(angle = -45,
                                   size = 7))+
  labs(x = "", y = "")+
  facet_rep_wrap(~study_period, repeat.tick.labels = TRUE)

各調査期間のCSIマトリックス間の相関は以下の通り。

# r <- matrix(ncol = 5, nrow = 5)
# 
# for(i in seq_along(sp)){
#   for(j in seq_along(sp)){
#     result <- mantel(CSI_sp_female[[i]], CSI_sp_female[[j]])
#     r[i,j] <- result$statistic
#   }
# }
# 
# rownames(r) <- sp
# colnames(r) <- sp
# 
# r

3 TY・ITとメスのCSIの計算

3.1 CSIの定義

ここでは、TY・ITとのCSIを算出する。

算出方法は Yamaguchi & Kazahari (2022) に従う。CSIの算出式は以下の通り。なお、\(i = 1,2,3,...,N\)はメスのIDを表す。\(P_i\)の算出の際にはTY/ITと毛づくろいしているポイントは除いた。対象としたのは2018年時点で5歳以上だった非発情の個体である。

  • \(G_{i}\): メスiとTY/ITの毛づくろい時間割合 (TYと毛づくろいしたポイント数/総ポイント数)
  • \(P_{i}\): メスiとTY/ITの近接時間割合 (3m以内にTYがいたポイント数/地上採食または休息したポイント数)
  • \(\overline{G}: \frac{1}{N} \sum^{N}_{i=1}G_{i}\)
  • \(\overline{P}: \frac{1}{N} \sum^{N}_{i=1}P_{i}\)

\[ CSI_{i} = \frac{\left(\frac{G_{i}}{\overline{G}} +\frac{P_{i}}{\overline{P}}\right)}{2} \;\;(i = 1,2,\dots,N) \]

使用するのは2018年交尾期、2019年交尾期、2019年非交尾期、2021年非交尾期、2022年非交尾期のデータである。

3.2 TYとメスのCSIの算出

3.2.1 データの読み込みと加工

まずは、TYの出入りがあった2019年交尾期2021年非交尾期についてTYが確認できた時間帯のデータを読み込む。

TY_19m <- read_excel("../Data/data/2019mating/2019mating_raw.xlsx",
                      sheet = "male_presence_long") %>% 
                filter(maleID == "TY") %>%
                mutate_at(5:6, ~make_datetime(year(date), month(date), mday(date),
                                              hour(.),minute(.))) %>% 
                mutate(date = as_date(date)) %>% 
                ## 分派集団にいるときは除く
                filter(groupID != 41 & groupID != 53 & groupID != 54)   

TY_21nm <- read_excel("../Data/data/2021nonmating/2021nonmating_raw.xlsx",
                      sheet = "male_presence_long") %>% 
                filter(maleID == "TY") %>%
                mutate_at(5:6, ~make_datetime(year(date), month(date), mday(date),
                                              hour(.),minute(.))) %>% 
                mutate(date = as_date(date))   

上記2つの期間の個体追跡データに、各個体追跡セッション中にTYが確認できていたかを表す列を追加し、TYが確認できたものだけを抽出する(focal_19m_bfocal_21nm_b)。

focal_19m %>% 
  distinct(no_focal, date, start_time, fin_time) %>% 
  left_join(TY_19m %>% select(date, male_presence, first, last), by = "date") %>% 
  ## 個体追跡セッション中にTYが確認されていたか  
  mutate(TY = ifelse((first <= start_time & last >= fin_time)|
                       (first >= start_time & first <= fin_time)|
                       (last >= start_time & last <= fin_time), 1, 0)) %>% 
  replace_na(list(TY = 0)) -> focal_19m_TY

focal_19m %>% 
  left_join(focal_19m_TY %>% select(no_focal, TY), by = "no_focal") %>% 
  filter(TY == "1") -> focal_19m_TY_b

focal_21nm %>% 
  distinct(no_focal, date, start_time, fin_time) %>% 
  left_join(TY_21nm %>% select(date, male_presence, first, last), by = "date") %>% 
  ## 個体追跡セッション中にTYが確認されていたか  
  mutate(TY = ifelse((first <= start_time & last >= fin_time)|
                       (first >= start_time & first <= fin_time)|
                       (last >= start_time & last <= fin_time), 1, 0)) %>% 
  replace_na(list(TY = 0)) -> focal_21nm_TY

focal_21nm %>% 
  left_join(focal_21nm_TY %>% select(no_focal, TY), by = "no_focal") %>% 
  filter(TY == "1") -> focal_21nm_TY_b

交尾期のデータについて、追跡個体が発情しているかの列を追加し、発情していない日のデータのみを抽出する(focal_18m_TY_bfocal_19m_TY_b)。

focal_18m %>% 
  left_join(female_18m %>% 
              rename(subject = femaleID) %>% 
              select(date, subject, rs2) , by = c("date", "subject")) %>% 
  filter(rs2 == "0") -> focal_18m_TY_b

focal_19m_TY_b %>% 
  left_join(female_19m %>% 
              rename(subject = femaleID) %>% 
              select(date, subject, rs2) , by = c("date", "subject")) %>% 
  filter(rs2 == "0") -> focal_19m_TY_b

3.2.2 CSIの算出(全調査期間)

3.2.2.1 毛づくろい時間割合の算出

算出のためにデータを加工する。

## 全期間を結合
focal_all_TY <- bind_rows(focal_18m_TY_b, focal_19m_TY_b, focal_19nm,
                          focal_21nm_TY_b, focal_22nm)

## 毛づくろい相手を表す列を追加  
focal_all_TY %>% 
  mutate(groom = ifelse(groomer == subject,groomee,groomer), 
         groom2 = ifelse(groomer2 == subject,groomee2,groomer2)) -> focal_all_TY_b

各メスの追跡時間(総瞬間サンプリングポイント数)を算出する。

focal_all_TY_b %>% 
  group_by(subject) %>% 
  summarise(dur = n()) %>% 
  ungroup() -> focal_duration_TY

各個体のTYとの毛づくろい時間を算出する。

focal_all_TY_b %>% 
  filter(groom == "TY"|groom2 == "TY") %>% 
  group_by(subject) %>% 
  summarise(no_groom = n()) %>% 
  ungroup()-> no_groom_TY

最後に、毛づくろい時間割合を算出する。

left_join(focal_duration_TY, no_groom_TY, by = "subject") %>% 
  replace_na(list(no_groom = 0)) %>% 
  mutate(prop_groom = no_groom/dur) -> prop_groom_TY

3.2.2.2 近接時間割合

地上採食/地上休息・毛づくろい時にTYが3m以内にいた時間割合を算出する。

まずは、地上採食/地上休息・毛づくろいのデータのみを抽出する。また、TYと毛づくろいしているポイントは除外する。

focal_all_TY_b %>% 
  filter(activity %in% c("F","R","G") & T_G == "G") %>% 
  select(no_focal, subject, x0_1m:x1_3m, groom, groom2, study_period) %>% 
  replace_na(list(groom = "NA", groom2 = "NA")) %>% 
  filter(!(groom == "TY"|groom2 == "TY")) -> focal_prox_TY

分母を算出する。なお、TYと毛づくろいしていたポイントは分母に含まないものとする。

focal_prox_TY %>% 
  group_by(subject) %>% 
  summarise(dur = n()) %>% 
  ungroup() -> prox_duration_TY

TYと3m以内に近接していたポイント数を算出する。

focal_prox_TY %>% 
  filter(str_detect(x0_1m, "TY")|str_detect(x1_3m, "TY")) %>% 
  group_by(subject) %>% 
  summarise(no_prox = n()) %>% 
  ungroup() -> no_prox_TY

最後に、近接時間割合を算出する。

left_join(prox_duration_TY, no_prox_TY, by = "subject") %>% 
  mutate(prop_prox = no_prox/dur) -> prop_prox_TY

3.2.2.3 CSIの算出

それでは、上記のデータを用いてCSIを算出する。

## 平均毛づくろい時間割合  
mean_groom_TY <-  mean(prop_groom_TY$prop_groom)

## 平均近接時間割合  
mean_prox_TY <- mean(prop_prox_TY$prop_prox)

## CSIの算出  
### 平均値  
prop_groom_TY %>% 
  left_join(prop_prox_TY, by = "subject") %>% 
  mutate(CSI_TY = ((prop_groom/mean_groom_TY) + (prop_prox/mean_prox_TY))/2) -> CSI_TY

算出した値が以下の通り。Kitは一日TYに長時間毛づくろいした日があり、それが強く影響したようだ。

CSI_TY %>% 
  mutate(subject = fct_reorder(subject, CSI_TY)) %>% 
  ggplot(aes(x = subject, y = CSI_TY))+
  geom_col(color = "black", fill = "white", size =0.3)+
  theme_classic(base_size=12, base_family = "Arial")+
  xlab("")+
  ylab(expression(paste(CSI, " with ", italic(TY))))+
  labs(title = "A")+
  coord_flip(ylim = c(0,5.3)) +
  scale_y_continuous(breaks = seq(0,5.2, by= 0.5),
                     expand= c(0,0))+
  theme_bw()+
  theme(aspect.ratio=1,
        axis.text.y = element_text(face = "bold.italic",
                                   hjust=0,
                                   family = "Times New Roman"),
        axis.text.x = element_text(family = "Times New Roman",
                                   face = "bold"),
        axis.ticks.y = element_blank(),
        axis.title.y = element_text(size = 8, 
                                    family = "Times New Roman"),
        axis.title.x = element_text(size = 10.5,
                                    family = "Times New Roman"),
        plot.margin=grid::unit(c(0,0,0,0), "mm"),
        plot.title = element_text(family = "Times New Roman",
                             hjust = 0)) -> p_CSI_TY

p_CSI_TY

# ggsave("figures/p_CSI_TY.tiff",
#        p_CSI_TY, width = 100, height = 100, dpi = 1500,
#        units = "mm")

3.2.3 CSIの算出(2018交尾期~2019交尾期)

3.2.3.1 毛づくろい時間割合の算出

算出のためにデータを加工する。

## 期間を結合
focal_fh_TY <- bind_rows(focal_18m_TY_b, focal_19m_TY_b, focal_19nm)

## 毛づくろい相手を表す列を追加  
focal_fh_TY %>% 
  mutate(groom = ifelse(groomer == subject,groomee,groomer), 
         groom2 = ifelse(groomer2 == subject,groomee2,groomer2)) -> focal_fh_TY_b

各メスの追跡時間(総瞬間サンプリングポイント数)を算出する。

focal_fh_TY_b %>% 
  group_by(subject) %>% 
  summarise(dur = n()) %>% 
  ungroup() -> focal_duration_fh_TY

各個体のTYとの毛づくろい時間を算出する。

focal_fh_TY_b %>% 
  filter(groom == "TY"|groom2 == "TY") %>% 
  group_by(subject) %>% 
  summarise(no_groom = n()) %>% 
  ungroup() -> no_groom_fh_TY

最後に、毛づくろい時間割合を算出する。

left_join(focal_duration_fh_TY, no_groom_fh_TY, by = "subject") %>% 
  replace_na(list(no_groom = 0)) %>% 
  mutate(prop_groom = no_groom/dur) -> prop_groom_fh_TY

3.2.3.2 近接時間割合

地上採食/地上休息・毛づくろい時にTYが3m以内にいた時間割合を算出する。

まずは、地上採食/地上休息・毛づくろいのデータのみを抽出する。また、TYと毛づくろいしているポイントは除外する。

focal_fh_TY_b %>% 
  filter(activity %in% c("F","R","G") & T_G == "G") %>% 
  select(no_focal, subject, x0_1m:x1_3m, groom, groom2) %>% 
  replace_na(list(groom = "NA", groom2 = "NA")) %>% 
  filter(!(groom == "TY"|groom2 == "TY")) -> focal_prox_fh_TY

分母を算出する。なお、TYと毛づくろいしていたポイントは分母に含まないものとする。

focal_prox_fh_TY %>% 
  group_by(subject) %>% 
  summarise(dur = n()) %>% 
  ungroup() -> prox_duration_fh_TY

TYと3m以内に近接していたポイント数を算出する。

focal_prox_fh_TY %>% 
  filter(str_detect(x0_1m, "TY")|str_detect(x1_3m, "TY")) %>% 
  group_by(subject) %>% 
  summarise(no_prox = n()) %>% 
  ungroup() -> no_prox_fh_TY

最後に、近接時間割合を算出する。

left_join(prox_duration_fh_TY, no_prox_fh_TY, by = "subject") %>% 
  mutate(prop_prox = no_prox/dur) -> prop_prox_fh_TY

3.2.3.3 CSIの算出

それでは、上記のデータを用いてCSIを算出する。

## 平均毛づくろい時間割合  
mean_groom_fh_TY <-  mean(prop_groom_fh_TY$prop_groom)

## 平均近接時間割合  
mean_prox_fh_TY <- mean(prop_prox_fh_TY$prop_prox)

## CSIの算出  
prop_groom_fh_TY %>% 
  left_join(prop_prox_fh_TY, by = "subject") %>% 
  mutate(CSI_TY = ((prop_groom/mean_groom_fh_TY) + (prop_prox/mean_prox_fh_TY))/2) -> CSI_fh_TY

算出した値が以下の通り。

CSI_fh_TY %>% 
  mutate(subject = fct_reorder(subject, CSI_TY)) %>% 
  ggplot(aes(x = subject, y = CSI_TY))+
  geom_col(color = "black", fill = "white", size =0.3)+
  theme_classic(base_size=12, base_family = "Arial")+
  xlab("")+
  ylab(expression(CSI[i]))+
  coord_flip(ylim = c(0,5.3)) +
  scale_y_continuous(breaks = seq(0,5.2, by= 0.5),
                     expand= c(0,0))+
  theme(aspect.ratio=1,
        axis.text.y = element_text(face = "italic",
                                   hjust=0),
        axis.ticks.y = element_blank(),
        axis.title.y = element_text(size = 8),
        axis.title.x = element_text(size = 10.5),
        plot.margin=grid::unit(c(0,0,0,0), "mm"),
        plot.title = element_text(hjust = 0.5)) 

3.2.4 CSIの算出(非交尾期のみ)

3.2.4.1 毛づくろい時間割合の算出

算出のためにデータを加工する。

## 全期間を結合
focal_nm_TY <- bind_rows(focal_19nm, focal_21nm_TY_b, focal_22nm)

## 毛づくろい相手を表す列を追加  
focal_nm_TY %>% 
  mutate(groom = ifelse(groomer == subject,groomee,groomer), 
         groom2 = ifelse(groomer2 == subject,groomee2,groomer2)) -> focal_nm_TY_b

各メスの追跡時間(総瞬間サンプリングポイント数)を算出する。

focal_nm_TY_b %>% 
  group_by(subject) %>% 
  summarise(dur = n()) %>% 
  ungroup() -> focal_duration_nm_TY

各個体のTYとの毛づくろい時間を算出する。

focal_nm_TY_b %>% 
  filter(groom == "TY"|groom2 == "TY") %>% 
  group_by(subject) %>% 
  summarise(no_groom = n()) %>% 
  ungroup() -> no_groom_nm_TY

最後に、毛づくろい時間割合を算出する。

left_join(focal_duration_nm_TY, no_groom_nm_TY, by = "subject") %>% 
  replace_na(list(no_groom = 0)) %>% 
  mutate(prop_groom = no_groom/dur) -> prop_groom_nm_TY

3.2.4.2 近接時間割合

地上採食/地上休息・毛づくろい時にTYが3m以内にいた時間割合を算出する。

まずは、地上採食/地上休息・毛づくろいのデータのみを抽出する。また、TYと毛づくろいしているポイントは除外する。

focal_nm_TY_b %>% 
  filter(activity %in% c("F","R","G") & T_G == "G") %>% 
  select(no_focal, subject, x0_1m:x1_3m, groom, groom2) %>% 
  replace_na(list(groom = "NA", groom2 = "NA")) %>% 
  filter(!(groom == "TY"|groom2 == "TY")) -> focal_prox_nm_TY

分母を算出する。なお、TYと毛づくろいしていたポイントは分母に含まないものとする。

focal_prox_nm_TY %>% 
  group_by(subject) %>% 
  summarise(dur = n()) %>% 
  ungroup() -> prox_duration_nm_TY

TYと3m以内に近接していたポイント数を算出する。

focal_prox_nm_TY %>% 
  filter(str_detect(x0_1m, "TY")|str_detect(x1_3m, "TY")) %>% 
  group_by(subject) %>% 
  summarise(no_prox = n()) %>% 
  ungroup() -> no_prox_nm_TY

最後に、近接時間割合を算出する。

left_join(prox_duration_nm_TY, no_prox_nm_TY, by = "subject") %>% 
  mutate(prop_prox = no_prox/dur) -> prop_prox_nm_TY

3.2.4.3 CSIの算出

それでは、上記のデータを用いてCSIを算出する。

## 平均毛づくろい時間割合  
mean_groom_nm_TY <-  mean(prop_groom_nm_TY$prop_groom)

## 平均近接時間割合  
mean_prox_nm_TY <- mean(prop_prox_nm_TY$prop_prox)

## CSIの算出  
prop_groom_nm_TY %>% 
  left_join(prop_prox_nm_TY, by = "subject") %>% 
  mutate(CSI_TY = ((prop_groom/mean_groom_nm_TY) + (prop_prox/mean_prox_nm_TY))/2) -> CSI_nm_TY

算出した値が以下の通り。

CSI_nm_TY %>% 
  mutate(subject = fct_reorder(subject, CSI_TY)) %>% 
  ggplot(aes(x = subject, y = CSI_TY))+
  geom_col(color = "black", fill = "white", size =0.3)+
  theme_classic(base_size=12, base_family = "Arial")+
  xlab("")+
  ylab(expression(CSI[i]))+
  coord_flip(ylim = c(0,5.3)) +
  scale_y_continuous(breaks = seq(0,5.2, by= 0.5),
                     expand= c(0,0))+
  theme(aspect.ratio=1,
        axis.text.y = element_text(face = "italic",
                                   hjust=0),
        axis.ticks.y = element_blank(),
        axis.title.y = element_text(size = 8),
        axis.title.x = element_text(size = 10.5),
        plot.margin=grid::unit(c(0,0,0,0), "mm"),
        plot.title = element_text(hjust = 0.5)) 

3.2.5 CSIの算出(調査期間ごと)

3.2.5.1 毛づくろい時間割合の算出

各メスの追跡時間(総瞬間サンプリングポイント数)を調査期間ごとに算出する。

focal_all_TY_b %>% 
  group_by(subject, study_period) %>% 
  summarise(dur = n()) %>% 
  ungroup() -> focal_duration_sp_TY

各個体のTYとの毛づくろい時間を算出する。

focal_all_TY_b %>% 
  filter(groom == "TY"|groom2 == "TY") %>% 
  group_by(subject, study_period) %>% 
  summarise(no_groom = n()) %>% 
  ungroup() -> no_groom_sp_TY

最後に、毛づくろい時間割合を算出する。

left_join(focal_duration_sp_TY, no_groom_sp_TY, by = c("subject","study_period")) %>% 
  replace_na(list(no_groom = 0)) %>% 
  mutate(prop_groom = no_groom/dur) -> prop_groom_sp_TY

3.2.5.2 近接時間割合

地上採食/地上休息・毛づくろい時にTYが3m以内にいた時間割合を算出する。

調査期間ごとに分母を算出する。なお、TYと毛づくろいしていたポイントは分母に含まないものとする。

focal_prox_TY %>% 
  group_by(subject, study_period) %>% 
  summarise(dur = n()) %>% 
  ungroup() -> prox_duration_sp_TY

TYと3m以内に近接していたポイント数を算出する。

focal_prox_TY %>% 
  filter(str_detect(x0_1m, "TY")|str_detect(x1_3m, "TY")) %>% 
  group_by(subject, study_period) %>%
  summarise(no_prox = n()) %>% 
  ungroup() -> no_prox_sp_TY

最後に、近接時間割合を算出する。

left_join(prox_duration_sp_TY, no_prox_sp_TY, by = c("subject","study_period")) %>% 
  replace_na(list(no_prox = 0)) %>% 
  mutate(prop_prox = no_prox/dur) -> prop_prox_sp_TY

3.2.5.3 CSIの算出

それでは、上記のデータを用いてCSIを算出する。

## CSIの算出  
prop_groom_sp_TY %>% 
  left_join(prop_prox_sp_TY, by = c("subject","study_period")) %>% 
  group_by(study_period) %>% 
  mutate(mean_groom = mean(prop_groom),
         mean_prox = mean(prop_prox)) %>% 
  ungroup() %>% 
  mutate(CSI_TY = ((prop_groom/mean_groom) + (prop_prox/mean_prox))/2) -> CSI_sp_TY

算出した値が以下の通り。

CSI_sp_TY %>% 
  ggplot(aes(x = study_period, y = CSI_TY))+
  geom_point()+
  geom_line(aes(group = subject))+
  theme_bw()+
  facet_rep_wrap(~subject, repeat.tick.labels = TRUE)+
  ylab(expression(CSI[i]))

3.3 ITとメスのCSIの算出

3.3.1 データの読み込みと加工

まずは、ITの出入りがあった2019年交尾期2021年非交尾期についてTYが確認できた時間帯のデータを読み込む。

IT_19m <- read_excel("../Data/data/2019mating/2019mating_raw.xlsx",
                      sheet = "male_presence_long") %>% 
                filter(maleID == "IT") %>%
                mutate_at(5:6, ~make_datetime(year(date), month(date), mday(date),
                                              hour(.),minute(.))) %>% 
                mutate(date = as_date(date)) %>% 
                ## 分派集団にいるときは除く
                filter(groupID != 41 & groupID != 53 & groupID != 54)   

IT_21nm <- read_excel("../Data/data/2021nonmating/2021nonmating_raw.xlsx",
                      sheet = "male_presence_long") %>% 
                filter(maleID == "IT") %>%
                mutate_at(5:6, ~make_datetime(year(date), month(date), mday(date),
                                              hour(.),minute(.))) %>% 
                mutate(date = as_date(date))   

上記2つの期間の個体追跡データに、各個体追跡セッション中にITが確認できていたかを表す列を追加し、ITが確認できたものだけを抽出する(focal_19m_bfocal_21nm_b)。

focal_19m %>% 
  distinct(no_focal, date, start_time, fin_time) %>% 
  left_join(IT_19m %>% select(date, male_presence, first, last), by = "date") %>% 
  ## 個体追跡セッション中にITが確認されていたか  
  mutate(IT = ifelse((first <= start_time & last >= fin_time)|
                       (first >= start_time & first <= fin_time)|
                       (last >= start_time & last <= fin_time), 1, 0)) %>% 
  replace_na(list(IT = 0)) -> focal_19m_IT

focal_19m %>% 
  left_join(focal_19m_IT %>% select(no_focal, IT), by = "no_focal") %>% 
  filter(IT == "1") -> focal_19m_IT_b

focal_21nm %>% 
  distinct(no_focal, date, start_time, fin_time) %>% 
  left_join(IT_21nm %>% select(date, male_presence, first, last), by = "date") %>% 
  ## 個体追跡セッション中にITが確認されていたか  
  mutate(IT = ifelse((first <= start_time & last >= fin_time)|
                       (first >= start_time & first <= fin_time)|
                       (last >= start_time & last <= fin_time), 1, 0)) %>% 
  replace_na(list(IT = 0)) -> focal_21nm_IT

focal_21nm %>% 
  left_join(focal_21nm_IT %>% select(no_focal, IT), by = "no_focal") %>% 
  filter(IT == "1") -> focal_21nm_IT_b

交尾期のデータについて、追跡個体が発情しているかの列を追加し、発情していない日のデータのみを抽出する(focal_18m_IT_bfocal_19m_IT_b)。

focal_18m %>% 
  left_join(female_18m %>% 
              rename(subject = femaleID) %>% 
              select(date, subject, rs2) , by = c("date", "subject")) %>% 
  filter(rs2 == "0") -> focal_18m_IT_b

focal_19m_IT_b %>% 
  left_join(female_19m %>% 
              rename(subject = femaleID) %>% 
              select(date, subject, rs2) , by = c("date", "subject")) %>% 
  filter(rs2 == "0") -> focal_19m_IT_b

3.3.2 CSIの算出(全調査期間)

3.3.2.1 毛づくろい時間割合の算出

算出のためにデータを加工する。

## 全期間を結合
focal_all_IT <- bind_rows(focal_18m_IT_b, focal_19m_IT_b, focal_19nm,
                          focal_21nm_IT_b)

## 毛づくろい相手を表す列を追加  
focal_all_IT %>% 
  mutate(groom = ifelse(groomer == subject,groomee,groomer), 
         groom2 = ifelse(groomer2 == subject,groomee2,groomer2)) -> focal_all_IT_b

各メスの追跡時間(総瞬間サンプリングポイント数)を算出する。

focal_all_IT_b %>% 
  group_by(subject) %>% 
  summarise(dur = n()) %>% 
  ungroup() -> focal_duration_IT

各個体のITとの毛づくろい時間を算出する。

focal_all_IT_b %>% 
  filter(groom == "IT"|groom2 == "IT") %>% 
  group_by(subject) %>% 
  summarise(no_groom = n())%>% 
  ungroup() -> no_groom_IT

最後に、毛づくろい時間割合を算出する。

left_join(focal_duration_IT, no_groom_IT, by = "subject") %>% 
  replace_na(list(no_groom = 0)) %>% 
  mutate(prop_groom = no_groom/dur) -> prop_groom_IT

3.3.2.2 近接時間割合

地上採食/地上休息・毛づくろい時にITが3m以内にいた時間割合を算出する。

まずは、地上採食/地上休息・毛づくろいのデータのみを抽出する。また、ITと毛づくろいしているポイントは除外する。

focal_all_IT_b %>% 
  filter(activity %in% c("F","R","G") & T_G == "G") %>% 
  select(no_focal, subject, x0_1m:x1_3m, groom, groom2, study_period) %>% 
  replace_na(list(groom = "NA", groom2 = "NA")) %>% 
  filter(!(groom == "IT"|groom2 == "IT")) -> focal_prox_IT

分母を算出する。なお、ITと毛づくろいしていたポイントは分母に含まないものとする。

focal_prox_IT %>% 
  group_by(subject) %>% 
  summarise(dur = n())%>% 
  ungroup() -> prox_duration_IT

ITと3m以内に近接していたポイント数を算出する。

focal_prox_IT %>% 
  filter(str_detect(x0_1m, "IT")|str_detect(x1_3m, "IT")) %>% 
  group_by(subject) %>% 
  summarise(no_prox = n())%>% 
  ungroup() -> no_prox_IT

最後に、近接時間割合を算出する。

left_join(prox_duration_IT, no_prox_IT, by = "subject") %>% 
  mutate(prop_prox = no_prox/dur) -> prop_prox_IT

3.3.2.3 CSIの算出

それでは、上記のデータを用いてCSIを算出する。

## 平均毛づくろい時間割合  
mean_groom_IT <-  mean(prop_groom_IT$prop_groom)

## 平均近接時間割合  
mean_prox_IT <- mean(prop_prox_IT$prop_prox)

## CSIの算出  
prop_groom_IT %>% 
  left_join(prop_prox_IT, by = "subject") %>% 
  mutate(CSI_IT = ((prop_groom/mean_groom_IT) + (prop_prox/mean_prox_IT))/2) -> CSI_IT

算出した値が以下の通り。Kitは一日ITに長時間毛づくろいした日があり、それが強く影響したようだ。

CSI_IT %>% 
  mutate(subject = fct_reorder(subject, CSI_IT)) %>% 
  ggplot(aes(x = subject, y = CSI_IT))+
  geom_col(color = "black", fill = "white", size =0.3)+
  theme_classic(base_size=12, base_family = "Times New Roman")+
  xlab("")+
  ylab(expression(paste(CSI, " with ", italic(IT))))+
  labs(title = "B")+
  coord_flip(ylim = c(0,6.0)) +
  scale_y_continuous(breaks = seq(0,7, by= 0.5),
                     expand= c(0,0))+
  theme_bw()+
  theme(aspect.ratio=1,
        axis.text.y = element_text(face = "bold.italic",
                                   hjust=0,
                                   family = "Times New Roman"),
        axis.text.x = element_text(family = "Times New Roman",
                                   face = "bold"),
        axis.ticks.y = element_blank(),
        axis.title.y = element_text(size = 8, 
                                    family = "Times New Roman"),
        axis.title.x = element_text(size = 10.5,
                                    family = "Times New Roman"),
        plot.margin=grid::unit(c(0,0,0,0), "mm"),
        plot.title = element_text(family = "Times New Roman",
                             hjust = 0)) -> p_CSI_IT

p_CSI_IT

# ggsave("figures/p_CSI_IT.tiff",
#        p_CSI_IT, width = 100, height = 100, dpi = 1500,
#        units = "mm")

3.3.3 CSIの算出(2018交尾期~2019交尾期)

3.3.3.1 毛づくろい時間割合の算出

算出のためにデータを加工する。

## 期間を結合
focal_fh_IT <- bind_rows(focal_18m_IT_b, focal_19m_IT_b, focal_19nm)

## 毛づくろい相手を表す列を追加  
focal_fh_IT %>% 
  mutate(groom = ifelse(groomer == subject,groomee,groomer), 
         groom2 = ifelse(groomer2 == subject,groomee2,groomer2)) -> focal_fh_IT_b

各メスの追跡時間(総瞬間サンプリングポイント数)を算出する。

focal_fh_IT_b %>% 
  group_by(subject) %>% 
  summarise(dur = n())%>% 
  ungroup() -> focal_duration_fh_IT

各個体のITとの毛づくろい時間を算出する。

focal_fh_IT_b %>% 
  filter(groom == "IT"|groom2 == "IT") %>% 
  group_by(subject) %>% 
  summarise(no_groom = n())%>% 
  ungroup() -> no_groom_fh_IT

最後に、毛づくろい時間割合を算出する。

left_join(focal_duration_fh_IT, no_groom_fh_IT, by = "subject") %>% 
  replace_na(list(no_groom = 0)) %>% 
  mutate(prop_groom = no_groom/dur) -> prop_groom_fh_IT

3.3.3.2 近接時間割合

地上採食/地上休息・毛づくろい時にITが3m以内にいた時間割合を算出する。

まずは、地上採食/地上休息・毛づくろいのデータのみを抽出する。また、ITと毛づくろいしているポイントは除外する。

focal_fh_IT_b %>% 
  filter(activity %in% c("F","R","G") & T_G == "G") %>% 
  select(no_focal, subject, x0_1m:x1_3m, groom, groom2) %>% 
  replace_na(list(groom = "NA", groom2 = "NA")) %>% 
  filter(!(groom == "IT"|groom2 == "IT")) -> focal_prox_fh_IT

分母を算出する。なお、ITと毛づくろいしていたポイントは分母に含まないものとする。

focal_prox_fh_IT %>% 
  group_by(subject) %>% 
  summarise(dur = n())%>% 
  ungroup() -> prox_duration_fh_IT

ITと3m以内に近接していたポイント数を算出する。

focal_prox_fh_IT %>% 
  filter(str_detect(x0_1m, "IT")|str_detect(x1_3m, "IT")) %>% 
  group_by(subject) %>% 
  summarise(no_prox = n())%>% 
  ungroup() -> no_prox_fh_IT

最後に、近接時間割合を算出する。

left_join(prox_duration_fh_IT, no_prox_fh_IT, by = "subject") %>% 
  replace_na(list(no_prox = 0)) %>% 
  mutate(prop_prox = no_prox/dur) -> prop_prox_fh_IT

3.3.3.3 CSIの算出

それでは、上記のデータを用いてCSIを算出する。

## 平均毛づくろい時間割合  
mean_groom_fh_IT <-  mean(prop_groom_fh_IT$prop_groom)

## 平均近接時間割合  
mean_prox_fh_IT <- mean(prop_prox_fh_IT$prop_prox)

## CSIの算出  
prop_groom_fh_IT %>% 
  left_join(prop_prox_fh_IT, by = "subject") %>% 
  mutate(CSI_IT = ((prop_groom/mean_groom_fh_IT) + (prop_prox/mean_prox_fh_IT))/2) -> CSI_fh_IT

算出した値が以下の通り。

CSI_fh_IT %>% 
  mutate(subject = fct_reorder(subject, CSI_IT)) %>% 
  ggplot(aes(x = subject, y = CSI_IT))+
  geom_col(color = "black", fill = "white", size =0.3)+
  theme_classic(base_size=12, base_family = "Arial")+
  xlab("")+
  ylab(expression(CSI[i]))+
  coord_flip(ylim = c(0,10)) +
  scale_y_continuous(breaks = seq(0,10, by= 1),
                     expand= c(0,0))+
  theme(aspect.ratio=1,
        axis.text.y = element_text(face = "italic",
                                   hjust=0),
        axis.ticks.y = element_blank(),
        axis.title.y = element_text(size = 8),
        axis.title.x = element_text(size = 10.5),
        plot.margin=grid::unit(c(0,0,0,0), "mm"),
        plot.title = element_text(hjust = 0.5)) 

3.3.4 CSIの算出(非交尾期のみ)

3.3.4.1 毛づくろい時間割合の算出

算出のためにデータを加工する。

## 全期間を結合
focal_nm_IT <- bind_rows(focal_19nm, focal_21nm_IT_b)

## 毛づくろい相手を表す列を追加  
focal_nm_IT %>% 
  mutate(groom = ifelse(groomer == subject,groomee,groomer), 
         groom2 = ifelse(groomer2 == subject,groomee2,groomer2)) -> focal_nm_IT_b

各メスの追跡時間(総瞬間サンプリングポイント数)を算出する。

focal_nm_IT_b %>% 
  group_by(subject) %>% 
  summarise(dur = n()) %>% 
  ungroup() -> focal_duration_nm_IT

各個体のITとの毛づくろい時間を算出する。

focal_nm_IT_b %>% 
  filter(groom == "IT"|groom2 == "IT") %>% 
  group_by(subject) %>% 
  summarise(no_groom = n()) %>% 
  ungroup() -> no_groom_nm_IT

最後に、毛づくろい時間割合を算出する。

left_join(focal_duration_nm_IT, no_groom_nm_IT, by = "subject") %>% 
  replace_na(list(no_groom = 0)) %>% 
  mutate(prop_groom = no_groom/dur) -> prop_groom_nm_IT

3.3.4.2 近接時間割合

地上採食/地上休息・毛づくろい時にITが3m以内にいた時間割合を算出する。

まずは、地上採食/地上休息・毛づくろいのデータのみを抽出する。また、ITと毛づくろいしているポイントは除外する。

focal_nm_IT_b %>% 
  filter(activity %in% c("F","R","G") & T_G == "G") %>% 
  select(no_focal, subject, x0_1m:x1_3m, groom, groom2) %>% 
  replace_na(list(groom = "NA", groom2 = "NA")) %>% 
  filter(!(groom == "IT"|groom2 == "IT")) -> focal_prox_nm_IT

分母を算出する。なお、ITと毛づくろいしていたポイントは分母に含まないものとする。

focal_prox_nm_IT %>% 
  group_by(subject) %>% 
  summarise(dur = n()) %>% 
  ungroup() -> prox_duration_nm_IT

ITと3m以内に近接していたポイント数を算出する。

focal_prox_nm_IT %>% 
  filter(str_detect(x0_1m, "IT")|str_detect(x1_3m, "IT")) %>% 
  group_by(subject) %>% 
  summarise(no_prox = n()) %>% 
  ungroup() -> no_prox_nm_IT

最後に、近接時間割合を算出する。

left_join(prox_duration_nm_IT, no_prox_nm_IT, by = "subject") %>% 
  replace_na(list(no_prox = 0)) %>% 
  mutate(prop_prox = no_prox/dur) -> prop_prox_nm_IT

3.3.4.3 CSIの算出

それでは、上記のデータを用いてCSIを算出する。

## 平均毛づくろい時間割合  
mean_groom_nm_IT <-  mean(prop_groom_nm_IT$prop_groom)

## 平均近接時間割合  
mean_prox_nm_IT <- mean(prop_prox_nm_IT$prop_prox)

## CSIの算出  
prop_groom_nm_IT %>% 
  left_join(prop_prox_nm_IT, by = "subject") %>% 
  mutate(CSI_IT = ((prop_groom/mean_groom_nm_IT) + (prop_prox/mean_prox_nm_IT))/2) -> CSI_nm_IT

算出した値が以下の通り。

CSI_nm_IT %>% 
  mutate(subject = fct_reorder(subject, CSI_IT)) %>% 
  ggplot(aes(x = subject, y = CSI_IT))+
  geom_col(color = "black", fill = "white", size =0.3)+
  theme_classic(base_size=12, base_family = "Arial")+
  xlab("")+
  ylab(expression(CSI[i]))+
  coord_flip(ylim = c(0,7)) +
  scale_y_continuous(breaks = seq(0,10, by= 1),
                     expand= c(0,0))+
  theme(aspect.ratio=1,
        axis.text.y = element_text(face = "italic",
                                   hjust=0),
        axis.ticks.y = element_blank(),
        axis.title.y = element_text(size = 8),
        axis.title.x = element_text(size = 10.5),
        plot.margin=grid::unit(c(0,0,0,0), "mm"),
        plot.title = element_text(hjust = 0.5)) 

3.3.5 CSIの算出(調査期間ごと)

3.3.5.1 毛づくろい時間割合の算出

各メスの追跡時間(総瞬間サンプリングポイント数)を調査期間ごとに算出する。

focal_all_IT_b %>% 
  group_by(subject, study_period) %>% 
  summarise(dur = n()) %>% 
  ungroup() -> focal_duration_sp_IT

各個体のITとの毛づくろい時間を算出する。

focal_all_IT_b %>% 
  filter(groom == "IT"|groom2 == "IT") %>% 
  group_by(subject, study_period) %>% 
  summarise(no_groom = n()) %>% 
  ungroup() -> no_groom_sp_IT

最後に、毛づくろい時間割合を算出する。

left_join(focal_duration_sp_IT, no_groom_sp_IT, by = c("subject","study_period")) %>% 
  replace_na(list(no_groom = 0)) %>% 
  mutate(prop_groom = no_groom/dur) -> prop_groom_sp_IT

3.3.5.2 近接時間割合

地上採食/地上休息・毛づくろい時にITが3m以内にいた時間割合を算出する。

調査期間ごとに分母を算出する。なお、ITと毛づくろいしていたポイントは分母に含まないものとする。

focal_prox_IT %>% 
  group_by(subject, study_period) %>% 
  summarise(dur = n()) %>% 
  ungroup() -> prox_duration_sp_IT

ITと3m以内に近接していたポイント数を算出する。

focal_prox_IT %>% 
  filter(str_detect(x0_1m, "IT")|str_detect(x1_3m, "IT")) %>% 
  group_by(subject, study_period) %>%
  summarise(no_prox = n()) %>% 
  ungroup() -> no_prox_sp_IT

最後に、近接時間割合を算出する。

left_join(prox_duration_sp_IT, no_prox_sp_IT, by = c("subject","study_period")) %>% 
  replace_na(list(no_prox = 0)) %>% 
  mutate(prop_prox = no_prox/dur) -> prop_prox_sp_IT

3.3.5.3 CSIの算出

それでは、上記のデータを用いてCSIを算出する。なお、平均毛づくろい時間割合が0の調査期間は分子の左側を1とした。

## CSIの算出  
prop_groom_sp_IT %>% 
  left_join(prop_prox_sp_IT, by = c("subject","study_period")) %>% 
  group_by(study_period) %>% 
  mutate(mean_groom = mean(prop_groom),
         mean_prox = mean(prop_prox)) %>% 
  ungroup() %>% 
  mutate(CSI_IT = ifelse(mean_groom != 0,
                         ((prop_groom/mean_groom) + (prop_prox/mean_prox))/2,
                         (1 + (prop_prox/mean_prox))/2)) -> CSI_sp_IT

算出した値が以下の通り。

CSI_sp_IT %>% 
  ggplot(aes(x = study_period, y = CSI_IT))+
  geom_point()+
  geom_line(aes(group = subject))+
  theme_bw()+
  facet_rep_wrap(~subject, repeat.tick.labels = TRUE)+
  ylab(expression(CSI[i]))

3.4 分母をTYとIT両方の頻度の平均値にする場合

CSIを算出する際の分母を、TYとITのデータの両方を合わせた平均にします。

3.4.1 CSIの算出

算出は以下のように行えます。

## 平均毛づくろい時間割合  
mean_groom_TYIT <-  mean(c(prop_groom_TY$prop_groom,
                           prop_groom_IT$prop_groom))

## 平均近接時間割合  
mean_prox_TYIT <- mean(c(prop_prox_TY$prop_prox,
                         prop_prox_IT$prop_prox))

## CSIの算出  
### 平均値  
prop_groom_TY %>% 
  left_join(prop_prox_TY, by = "subject") %>% 
  mutate(CSI_TY = ((prop_groom/mean_groom_TYIT) + (prop_prox/mean_prox_TYIT))/2) -> CSI_TY_combined

prop_groom_IT %>% 
  left_join(prop_prox_IT, by = "subject") %>% 
  mutate(CSI_IT = ((prop_groom/mean_groom_TYIT) + (prop_prox/mean_prox_TYIT))/2) -> CSI_IT_combined  

3.4.2 図示

図示すると以下のようになります。

CSI_TY_combined %>% 
  mutate(subject = fct_reorder(subject, CSI_TY)) %>% 
  ggplot(aes(x = subject, y = CSI_TY))+
  geom_col(color = "black", fill = "white", size =0.3)+
  theme_classic(base_size=12, base_family = "Arial")+
  xlab("")+
  ylab(expression(paste(CSI, " with ", italic(TY))))+
  labs(title = "(a)")+
  coord_flip(ylim = c(0,5.15)) +
  scale_y_continuous(breaks = seq(0,5.2, by= 0.5),
                     expand= c(0,0))+
  theme_bw()+
  theme(aspect.ratio=1,
        axis.text.y = element_text(face = "italic",
                                   hjust=0,
                                   family = "Times New Roman"),
        axis.text.x = element_text(family = "Times New Roman"),
        axis.ticks.y = element_blank(),
        axis.title.y = element_text(size = 8, 
                                    family = "Times New Roman"),
        axis.title.x = element_text(size = 10.5,
                                    family = "Times New Roman"),
        plot.margin=grid::unit(c(0,0,0,0), "mm"),
        plot.title = element_text(family = "Times New Roman",
                             hjust = 0)) -> p_CSI_TY_combined

# ggsave("figures/p_CSI_TY_combined.tif",
#        p_CSI_TY_combined, width = 100, height = 100, dpi = 2000,
#        units = "mm")

CSI_IT_combined %>% 
  mutate(subject = fct_reorder(subject, CSI_IT)) %>% 
  ggplot(aes(x = subject, y = CSI_IT))+
  geom_col(color = "black", fill = "white", size =0.3)+
  theme_classic(base_size=12, base_family = "Arial")+
  xlab("")+
  ylab(expression(paste(CSI, " with ", italic(IT))))+
  labs(title = "(b)")+
  coord_flip(ylim = c(0,4)) +
  scale_y_continuous(breaks = seq(0,4, by= 0.5),
                     expand= c(0,0))+
  theme_bw()+
  theme(aspect.ratio=1,
        axis.text.y = element_text(face = "italic",
                                   hjust=0,
                                   family = "Times New Roman"),
        axis.text.x = element_text(family = "Times New Roman"),
        axis.ticks.y = element_blank(),
        axis.title.y = element_text(size = 8, 
                                    family = "Times New Roman"),
        axis.title.x = element_text(size = 10.5,
                                    family = "Times New Roman"),
        plot.margin=grid::unit(c(0,0,0,0), "mm"),
        plot.title = element_text(family = "Times New Roman",
                             hjust = 0)) -> p_CSI_IT_combined

# ggsave("figures/p_CSI_IT_combined.tif",
#        p_CSI_IT_combined, width = 100, height = 100, dpi = 2000,
#        units = "mm")

p_CSI_TYIT_combined <- p_CSI_TY_combined + p_CSI_IT_combined

# ggsave("figures/p_CSI_TYIT_combined.tif", p_CSI_TYIT_combined,
#        width = 180, height = 105, units = "mm", dpi = 2000)

3.5 算出に使用したメスの個体追跡データの概要

以下に、全期間のCSIを算出するのに使用した各メスの個体追跡データの総時間(総瞬間サンプリングポイント数)の情報をまとめる。

focal_duration_TY %>% 
  left_join(focal_duration_IT, by = "subject") %>% 
  rename(femaleID = 1, "Total points(TY)" = 2, "Total points (IT)" = 3) %>% 
  flextable() %>% 
  colformat_double(digits=2) %>% 
  set_table_properties(layout="autofit",width = 1) %>% 
  autofit(add_w = 0.2) %>% 
  font(part = "header", fontname = "Times New Roman") %>% 
  font(part = "body", j=1:3, fontname = "Times New Roman") %>% 
  theme_zebra() %>% 
  hline_bottom() %>% 
  hline_top() %>% 
  align(j=2:3,part = "all",align = "center") -> table_focal_duration

table_focal_duration

femaleID

Total points(TY)

Total points (IT)

Aka

2,306

1,575

Ako

2,459

1,704

Hen

2,621

1,689

Hot

2,268

1,603

Kil

2,561

1,655

Kit

2,302

1,545

Koh

2,401

1,731

Kol

2,627

1,715

Kun

2,099

1,344

Kur

1,120

1,120

Mal

2,458

1,640

Mei

2,742

1,727

Mif

1,225

547

Mik

2,413

1,717

Ntr

2,800

1,801

Tam

1,458

1,458

Ten

2,618

1,726

Tot

2,720

1,848

4 個体追跡データの加工

以下では、個体追跡データを分析のために加工します。

4.1 必要なデータの読み込み

以下、必要なデータを読み込みます。

## メスを確認した時間  
female_time <- read_csv("../Data/data/others/female_pre_time.csv") %>% 
  mutate(date = as_date(date))

#個体の属性情報
att <- read_csv("../Data/data/others/attributes_sp_over6.csv")

4.2 データの加工

4.2.1 メスの情報

まず、個体追跡データに以下のような加工を行い、新たに列を作成します。

  • groupID: 追跡メスがどの集団にいたか
  • max_female: 最大メス数
  • no_female: 確認メス数
  • no_est: 発情メス数
## 各フォーカル個体がどの集団で観察されたか  
female_presence_group <- group_all %>% 
  pivot_longer(Kur:Yun,
               names_to = "subject",
               values_to = "presence") %>% 
  filter(presence == 1) %>% 
  ## 2020年11月13の2つ目の集団は除く
  filter(!groupID %in% c("m20_52")) %>% 
  select(date, groupID, subject) 

## 各観察日の最大個体  
max_female <- female_time %>% 
  pivot_longer(cols = Kur:Cur,
               names_to = "femaleID",
               values_to = "presence") %>% 
  filter(is.na(presence)|(presence != "DD" & presence != "NS")) %>% 
  group_by(date) %>% 
  summarise(max_female = n()) 

## 各観察日の確認メス数(6歳以上)と発情メス数
no_female <- group_all %>% 
  select(groupID, study_period, date, Kur:Yun) %>% 
  select(-c(TY,IT, LK, KR, KM, TG)) %>% 
  pivot_longer(cols = Kur:Yun,
               names_to = "femaleID",
               values_to = "presence") %>% 
  left_join(att) %>% 
  ## 6歳以上の個体のみを抽出
  filter(age >= 6) %>% 
  left_join(female_all %>% select(date, femaleID, rs2)) %>% 
  ## groupIDごとに個体数を算出
  group_by(date, groupID, study_period) %>% 
  summarise(no_female = sum(presence, na.rm = TRUE),
            no_est = sum(rs2, na.rm = TRUE)) %>% 
  ungroup()
  
## 結合  
focal_combined_b <- focal_combined %>% 
  mutate(date = as_date(date)) %>% 
  left_join(female_presence_group,
            by = c("subject", "date")) %>% 
  left_join(max_female,
            by = c("date")) %>% 
  left_join(no_female)

4.2.2 オスの情報

続いて、交尾期については追跡日に確認された群れ外オス数、オス数を追加します。

まず、生データを用いて群れ外オス数とオス数を算出します。

no_ntm_m18 <- male_18m %>% 
  mutate(ntm = if_else(maleID %in% c("TY", "IT", "KR", "LK"),
                       0, 1)) %>% 
  group_by(date) %>% 
  summarise(no_male = sum(presence, na.rm = TRUE),
            no_ntm = sum(ntm == 1 & presence == 1, na.rm = TRUE))

no_ntm_m19 <- male_19m %>% 
  mutate(ntm = if_else(maleID %in% c("TY", "IT", "KR", "LK"),
                       0, 1)) %>% 
  group_by(date) %>% 
  summarise(no_male = sum(presence, na.rm = TRUE),
            no_ntm = sum(ntm == 1 & presence == 1, na.rm = TRUE))

no_ntm_m20 <- male_20m %>% 
  mutate(ntm = if_else(maleID %in% c("TY", "IT", "KR", "LK", "KM"),
                       0, 1)) %>% 
  group_by(date) %>% 
  summarise(no_male = sum(presence, na.rm = TRUE),
            no_ntm = sum(ntm == 1 & presence == 1, na.rm = TRUE))

no_ntm_m21 <- male_21m %>% 
  mutate(ntm = if_else(maleID %in% c("TY", "IT", "KR", "LK", "KM", "TG"),
                       0, 1)) %>% 
  group_by(date) %>% 
  summarise(no_male = sum(presence, na.rm = TRUE),
            no_ntm = sum(ntm == 1 & presence == 1, na.rm = TRUE))

sum_ntm <- bind_rows(no_ntm_m18, no_ntm_m19, no_ntm_m20, no_ntm_m21)


これを結合します。

focal_combined_c <- focal_combined_b %>% 
  left_join(sum_ntm)

4.2.3 TYとITの確認状況

調査期間中、第一位オスのタイヨウ(TY)と第二位オスのイツモ(IT)の群れへの出入りが頻繁に観察されました。
そこで、彼らが個体追跡時にいたか否かについての列(TYIT)を追加します。

まず、彼らが確認できた時刻に関するデータを読み込む。

TYIT_presence_time <- read_excel("../Data/data/2019mating/2019mating_raw.xlsx",
                                 sheet = "male_presence_long") %>% 
  mutate(study_period = "m19") %>% 
  mutate(groupID = str_c(study_period,"_", groupID)) %>% 
  select(maleID, date, groupID, male_presence, first, last) %>% 
  bind_rows(read_excel("../Data/data/2020mating/2020mating_raw.xlsx",
                       sheet = "male_presence_long") %>% 
              mutate(study_period = "m20") %>% 
              mutate(groupID = str_c(study_period,"_", groupID)) %>% 
              select(maleID, date, groupID, male_presence, first, last)) %>% 
  bind_rows(read_excel("../Data/data/2021mating/2021mating_raw.xlsx",
                       sheet = "male_presence_long") %>% 
              mutate(study_period = "m21") %>% 
              mutate(groupID = str_c(study_period,"_", groupID)) %>% 
              select(maleID, date, groupID, male_presence, first, last)) %>% 
  bind_rows(read_excel("../Data/data/2019nonmating/2019nonmating_raw.xlsx",
                       sheet = "male_presence_long") %>% 
              mutate(study_period = "nm19") %>% 
              mutate(groupID = str_c(study_period,"_", groupID)) %>% 
              select(maleID, date, groupID, male_presence, first, last)) %>% 
  bind_rows(read_excel("../Data/data/2020nonmating/2020nonmating_raw.xlsx",
                       sheet = "male_presence_long") %>% 
              mutate(study_period = "nm20") %>% 
              mutate(groupID = str_c(study_period,"_", groupID)) %>% 
              select(maleID, date, groupID, male_presence, first, last)) %>% 
  bind_rows(read_excel("../Data/data/2021nonmating/2021nonmating_raw.xlsx",
                       sheet = "male_presence_long") %>% 
              mutate(study_period = "nm21") %>% 
              mutate(groupID = str_c(study_period,"_", groupID)) %>% 
              select(maleID, date, groupID, male_presence, first, last)) %>% 
  bind_rows(read_excel("../Data/data/2022nonmating/2022nonmating_raw.xlsx",
                       sheet = "male_presence_long") %>% 
              mutate(study_period = "nm22") %>% 
              mutate(groupID = str_c(study_period,"_", groupID)) %>% 
              select(maleID, date, groupID, male_presence, first, last)) %>% 
  mutate(across(c(first, last),
                ~make_datetime(year(date), month(date), mday(date), hour(.),minute(.)))) 

## 横長にする
TYIT_presence_time_wide <- TYIT_presence_time %>% 
  filter(maleID == "TY") %>% 
  rename(first_TY = first, 
         last_TY = last,
         TY = male_presence) %>% 
  select(-maleID) %>% 
  left_join(TYIT_presence_time %>% 
              filter(maleID == "IT") %>% 
              rename(first_IT = first, 
                     last_IT = last,
                     IT = male_presence) %>% 
              select(-maleID)) %>% 
  mutate(date = as_date(date))

TYIT_presence_time_wide


これを利用し、個体追跡中にTYまたはITがいたと考えられるかを表す列を作成します。

  • TY_presence/IT_presence
    • 1: 個体追跡開始前または個体追跡中にTYを確認、かつTYがいなくなったと思われる前に開始
    • 0: その他
focal_combined_c %>% 
  left_join(group_all %>% select(date, groupID, TY, IT)) %>% 
  left_join(TYIT_presence_time_wide %>% select(-TY, -IT)) %>% 
  mutate(TY_presence = case_when(
    study_period == "m18" ~ TY,
    is.na(first_TY) ~ 0,
    start_time >= first_TY & start_time <= last_TY ~ 1,
    fin_time >= first_TY ~ 1,
    start_time >= last_TY ~ 0,
    fin_time <= first_TY ~ 0,
    .default = 1
  )) %>% 
  mutate(IT_presence = case_when(
    study_period == "m18" ~ IT,
    is.na(first_IT) ~ 0,
    start_time >= first_IT & start_time <= last_IT ~ 1,
    fin_time >= first_IT ~ 1,
    start_time >= last_IT ~ 0,
    fin_time <= first_IT ~ 0,
    .default = 1
  )) -> focal_combined_d

4.2.4 フォーカルを攻撃したオスのID

最後に、フォーカルが攻撃を受けた際の攻撃したオスのIDを表す列を作成します。

focal_combined_d %>% 
  ## NAになっているところを0に
  replace_na(list(agg_focal = 0)) %>% 
  replace_na(list(victim1 = "NA", victim2 = "NA", victim3 = "NA",
                  aggressor1 = "NA", aggressor2 = "NA", aggressor = "NA")) %>% 
  mutate(if_victim1 = ifelse(agg_focal == 1 & str_detect(victim1, subject), 1, 0),
         if_victim2 = ifelse(agg_focal == 1 & str_detect(victim2, subject), 1, 0),
         if_victim3 = ifelse(agg_focal == 1 & str_detect(victim3, subject), 1, 0)) %>% 
  mutate(aggressor_focal = case_when(
    if_victim1 == 1 & if_victim2 == 1 & if_victim3 == 1 ~ str_c(aggressor1,aggressor2,aggressor3,sep =","),
    if_victim1 == 1 & if_victim2 == 1 ~ str_c(aggressor1, aggressor2,sep =","),
    if_victim1 == 1 & if_victim3 == 1 ~ str_c(aggressor1, aggressor3,sep =","),
    if_victim2 == 1 & if_victim3 == 1 ~ str_c(aggressor2, aggressor3,sep =","),
    if_victim1 == 1 ~ aggressor1,
    if_victim2 == 1 ~ aggressor2,
    if_victim3 == 1 ~ aggressor3,
    .default = NA
  )) %>% 
  ## aggressor2には、追跡個体が同時刻に2頭から攻撃を受けたときのみ値が入る  
  separate(aggressor_focal, into = c("aggressor_focal1","aggressor_focal2"), 
           sep = ",") -> focal_combined_final

加工したデータは以下の通り。

datatable(focal_combined_final,
          options = list(scrollX = 60),
          filter = list(position ="top"))

4.2.5 フォーカルリスト

フォーカル個体のリストを作成します。

focal_list <- focal_combined_final %>% 
  mutate(no_focal = str_c(study_period, "_", no_focal)) %>% 
  group_by(no_focal, date, subject) %>% 
  summarise(n_points = n()) %>% 
  ungroup()

focal_list

5 オスの攻撃と怪我の関連

本章では、オスから受けた攻撃の頻度と怪我の有無の関連について検討します。

5.1 データの加工

まず怪我のデータを読み込みます。

## 怪我データの読み込み  
injury_m18 <- read_excel("../Data/data/2018mating/2018mating_raw.xlsx", sheet = "injury")
injury_m19 <- read_excel("../Data/data/2019mating/2019mating_raw.xlsx", sheet = "injury")
injury_m20 <- read_excel("../Data/data/2020mating/2020mating_raw.xlsx",sheet = "injury")
injury_m21 <- read_excel("../Data/data/2021mating/2021mating_raw.xlsx",sheet = "injury")
injury_nm19 <- read_excel("../Data/data/2019nonmating/2019nonmating_raw.xlsx", sheet = "injury")
injury_nm21 <- read_excel("../Data/data/2021nonmating/2021nonmating_raw.xlsx",sheet = "injury")
injury_nm22 <- read_excel("../Data/data/2022nonmating/2022nonmating_raw.xlsx",sheet = "injury")

injury_all <- bind_rows(injury_m18,injury_m19,injury_m20,injury_m21,
                        injury_nm19,injury_nm21,injury_nm22) %>% 
  pivot_longer(Kur:Har, names_to = "femaleID", values_to = "injury") %>% 
  mutate(injury01 = ifelse(injury >= 1, 1, 0)) %>% 
  mutate(date = as_date(date))

続いて、攻撃データを加工し、各個体の被攻撃頻度を日ごとに算出します。

### 個体ごとの被攻撃頻度  
agg_each <- aggression_all %>% 
  left_join(att, by = c("study_period", "femaleID")) %>% 
  filter(!is.na(age)) %>% 
  group_by(date, femaleID) %>% 
  summarise(no_agg = n()) %>% 
  ungroup() %>% 
  complete(date,femaleID) %>% 
  replace_na(list(no_agg = 0)) %>% 
  mutate(study_period = str_c("m",str_sub(year(date),3,4))) %>% 
  left_join(att, by = c("study_period", "femaleID")) %>% 
  ## 6歳以上を抽出  
  filter(!is.na(age)) 

### 個体データに結合  
## 各集団で確認できたメスと集団を追跡した時間を算出
group_all %>% 
  mutate(date = as_date(date),
         duration1 = as.numeric(difftime(fin, start, units = "mins")),
         duration2 = as.numeric(difftime(restart, suspend, units = "mins")))%>%
  mutate(duration1 = replace_na(duration1,0),
         duration2 = replace_na(duration2,0)) %>%
  mutate(duration = as.numeric(duration1 - duration2)) %>%  
  select(groupID, study_period, date, duration, Kur:Yun) %>% 
  select(-c(TY,IT, LK, KR, KM, TG)) %>% 
  pivot_longer(cols = Kur:Yun,
               names_to = "femaleID",
               values_to = "presence") %>% 
  left_join(att) %>% 
  filter(presence == 1 & age >= 6)  %>% 
  filter(date >= "2018-10-08") %>% 
  left_join(agg_each %>% select(date, femaleID, no_agg, study_period)) %>% 
  replace_na(list(no_agg = 0)) -> agg_rate_each

最後に、被攻撃頻度のデータと怪我データを結合します。また、フォーカルをその日にしたか否かの列も作成します。

## データを結合  
injury_final <- agg_rate_each %>% 
  left_join(injury_all, by = c("date", "femaleID")) %>% 
  left_join(female_all %>% select(date, femaleID, rs2),
            by = c("date", "femaleID")) %>%  
  drop_na(injury) %>% 
  left_join(focal_list %>% 
              group_by(date, subject) %>% 
              summarise(n_points = sum(n_points)) %>% 
              ungroup() %>% 
              select(date, subject, n_points) %>% 
              rename(femaleID = subject) %>% 
              mutate(focal = 1)) %>% 
  replace_na(list(n_points = 0,
                  focal = 0))

以下のようになりました。

datatable(injury_final)


5.2 分析

5.2.1 交尾期と非交尾期の比較

まず、交尾期と非交尾期の怪我頻度の違いを調べます。

5.2.1.1 モデリング

分析には、十分に観察を行った日(観察時間が240分以上)のみを用います。

  • 応答変数: 怪我が確認されたか否か(injury01)
  • 説明変数: 季節(mating_an vs mating_es vs nonmating)、観察時間(duration)
  • ランダム切片: メスID(femaleID)
injury_analysis_season <- injury_final %>% 
  filter(duration >= 240) %>% 
  mutate(season = if_else(str_detect(study_period, "nm"),
                          "nonmating", "mating")) %>% 
  mutate(season2 = case_when(
    season == "nonmating" ~ "nonmating",
    season == "mating" & rs2 == 1 ~ "mating_es",
    season == "mating" & rs2 == 0 ~ "mating_an",
    .default = NA
  ))

m_injury_season <- brm(injury01 ~ season2 + (1|femaleID),
                       family = bernoulli,
                       iter = 11000, warmup = 1000, seed = 13,
                       prior = c(prior(student_t(4,0,15), class = "b"),
                                 prior(student_t(4,0,15), class = "Intercept"),
                                 prior(student_t(4,0,10), class = "sd")),
                       control=list(adapt_delta = 0.9999, max_treedepth = 20),
                       backend = "cmdstanr",
                       data = injury_analysis_season,
                       file = "model/m_injury_season.rds")

5.2.1.2 モデルチェック

まず、DHARMaパッケージ(Hartig, 2022)とDHARMa.helperパッケージ(Rodríguez-Sánchez, 2023)でモデルの前提が満たされているかを確認する。いずれのモデルも特に問題はないよう。

## 全期間    
dh_injury_season <- dh_check_brms(m_injury_season, quantreg = TRUE)


ゼロ過剰などの問題もない。

testZeroInflation(dh_injury_season)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1, p-value = 1
## alternative hypothesis: two.sided


bayesplotパッケージ(Gabry & Mahr, 2022)pp_check関数で、事後分布からの予測分布と実測値の分布を比較しても大きな乖離はない。

## 全期間    
pp_check(m_injury_season, ndraws = 100)+
  theme_bw()+
  theme(aspect.ratio = 0.9)


多重共線性のチェックも行ったが、VIFに問題はない。

## 全期間  
check_collinearity(m_injury_season)
## NULL


Rhatにも問題はなく、収束の問題はないよう。

data.frame(Rhat = brms::rhat(m_injury_season)) %>% 
  ggplot(aes(x = Rhat))+
  geom_histogram(fill = "white",
                 color = "black")+
  theme_bw()+
  theme(aspect.ratio = 1)


5.2.1.3 結果の確認

5.2.1.3.1 推定された係数
まずは、推定された係数を確認します。


5.2.1.3.2 多重比較

続いて、季節/発情の有無について多重比較を行います。いずれも有意な差があることが分かります。

emm_injury_season <- emmeans(m_injury_season,
                             ~ season2,
                             type = "response") 

pairs(emm_injury_season)
##  contrast              odds.ratio lower.HPD upper.HPD
##  mating_an / mating_es      0.184     0.136      0.24
##  mating_an / nonmating      2.096     1.393      2.93
##  mating_es / nonmating     11.376     7.423     16.36
## 
## Point estimate displayed: median 
## Results are back-transformed from the log odds ratio scale 
## HPD interval probability: 0.95

5.2.1.4 結果の図示

最後に、結果を図示します。

injury_analysis_season %>% 
  group_by(femaleID, season2) %>% 
  summarise(prop_injury = mean(injury01)) %>% 
  mutate(season2 = fct_relevel(season2,
                               "nonmating", "mating_an")) %>% 
  ggplot(aes(x = season2, y = prop_injury)) +
  geom_violin(bw = 0.03,
              scale = "width")+
  geom_boxplot(size = 0.5,
               linewidth = 0.4,
               width = 0.2,
               fill = "grey",
               outlier.alpha = 0) +
  geom_point(data = emm_injury_season %>% 
               data.frame() %>% 
               mutate(season2 = fct_relevel(season2,
                                            "nonmating", "mating_an")),
             aes(y = response),
             size = 3,
             color = "white") +
  geom_errorbar(data = emm_injury_season %>% 
               data.frame() %>% 
               mutate(season2 = fct_relevel(season2,
                                            "nonmating", "mating_an")),
             aes(y = response, ymin = lower.HPD,
                 ymax = upper.HPD),
             size = 0.4,
             color = "white",
             width = 0.1) +
  geom_signif(comparisons = list(c("nonmating", "mating_es"),
                                 c("nonmating", "mating_an"),
                                 c("mating_an", "mating_es")),
              y_position = c(0.55, 0.1, 0.5),
              annotations = rep("***", 3)) +
  theme_bw(base_size = 12) +
  scale_x_discrete(labels = c("NMS", "MS (anestrus)", "MS (estrous)")) +
  labs(x = "", y = "Proportion of injured days (n/days)") +
  theme(text = element_text(family = "Times New Roman"),
        aspect.ratio = 1) -> p_injury_season

p_injury_season

ggsave("figures/p_injury_season.tif", p_injury_season,
       dpi = 2000, units = "mm",
       width = 100, height = 100)


5.2.2 被攻撃頻度と怪我の有無の関連

続いて、交尾期のデータを対象にメスの被攻撃頻度と怪我の有無の間に関連があるかを調べます。

5.2.2.1 モデリング

5.2.3 モデリング

それでは、モデリングを行います。分析には、その日個体追跡を行わなかった個体のデータのみを用います。これは、個体追跡をした個体は被攻撃頻度が他の個体に比べて過大評価されてしまう傾向があるからです。また、先ほどと同様に240分以上群れを追跡した日のデータのみを用います。

  • 分布: ベルヌーイ分布
  • リンク関数: ロジット関数
  • 応答変数: 怪我が確認されたか否か(injury01)
  • 説明変数: オスからの攻撃頻度(rate_agg)、メスの発情の有無(rs2)、追跡時間(duration)、 調査期間(study_period)
  • ランダム切片: メスID(femaleID)
injury_analysis <- injury_final %>% 
  filter(duration >= 240) %>% 
  filter(focal == 0) %>% 
  filter(!str_detect(study_period, "nm")) %>% 
  mutate(rate_agg = no_agg*60/duration) %>% 
  mutate(rate_agg_std = standardize(rate_agg))

m_agginj <- brm(injury01 ~ rate_agg_std*rs2  + study_period + (1|femaleID),
                family = bernoulli,
                iter = 11000, warmup = 1000, seed = 122,
                 prior = c(prior(student_t(4,0,15), class = "b"),
                           prior(student_t(4,0,15), class = "Intercept"),
                           prior(student_t(4,0,10), class = "sd")),
                 control=list(adapt_delta = 0.9999, max_treedepth = 20),
                 backend = "cmdstanr",
                 data = injury_analysis,
                 file = "model/m_agginj.rds")

5.2.3.1 モデルチェック

まず、DHARMaパッケージ(Hartig, 2022)とDHARMa.helperパッケージ(Rodríguez-Sánchez, 2023)でモデルの前提が満たされているかを確認する。いずれのモデルも特に問題はないよう。

## 全期間    
dh_agginj <- dh_check_brms(m_agginj, quantreg = TRUE)


ゼロ過剰などの問題もない。

testZeroInflation(dh_agginj)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0001, p-value = 1
## alternative hypothesis: two.sided


bayesplotパッケージ(Gabry & Mahr, 2022)pp_check関数で、事後分布からの予測分布と実測値の分布を比較しても大きな乖離はない。

## 全期間    
pp_check(m_agginj, ndraws = 100)+
  theme_bw()+
  theme(aspect.ratio = 0.9)


多重共線性のチェックも行ったが、VIFに問題はない。

## 全期間  
check_collinearity(m_agginj)


Rhatにも問題はなく、収束の問題はないよう。

data.frame(Rhat = brms::rhat(m_agginj)) %>% 
  ggplot(aes(x = Rhat))+
  geom_histogram(fill = "white",
                 color = "black")+
  theme_bw()+
  theme(aspect.ratio = 1)


5.2.3.2 結果の確認

5.2.3.2.1 推定された係数
まずは、推定された係数を確認します。


5.2.3.2.2 交互作用項の検討

続いて、季節/発情の有無について多重比較を行います。いずれも有意な差があることが分かります。

estimate_slopes(m_agginj,
                trend = "rate_agg_std",
                by = "rs2 = c(0,1)",
                ci_method = "eti") %>% 
  data.frame()

5.2.3.3 結果の図示

predict_agginj <- predict_response(m_agginj,
                                   terms = c("rate_agg_std[-0.4:13.4, by = 0.1]", "rs2[0:1, by = 1]"),
                                   type = "fixed",
                                   margin = "marginalmeans",
                                   interval = "confidence",
                                   ci_method = "ci",
                                   ci_level = 0.95) %>% 
  data.frame() %>% 
  rename(rate_agg_std = x,
         rs2 = group) %>% 
  mutate(rate_agg = rate_agg_std*sd(injury_analysis$rate_agg) + mean(injury_analysis$rate_agg)) %>% 
  mutate(rs2 = as.factor(rs2))

injury_analysis %>% 
  mutate(rs2 = as.factor(rs2)) %>% 
  ggplot(aes(x = rate_agg, y = injury01))+
  geom_point(shape = "|", size = 3,
             alpha = 0.5,
             aes(color = rs2)) +
  scale_color_nejm(labels = c("Anestrus", "Estrous")) +
  geom_line(data = predict_agginj,
            aes(y = predicted,
                linetype = rs2),
            linewidth = 0.4) +
  geom_ribbon(data = predict_agginj,
              aes(y = predicted,
                  ymax = conf.high, ymin = conf.low,
                  fill = rs2),
              alpha = 0.2) +
  scale_fill_nejm(labels = c("Anestrus", "Estrous")) +
  scale_linetype_manual(values = c("solid", "dashed"),
                        labels = c("Anestrus", "Estrous")) +
  theme_bw(base_size = 12) +
  scale_y_continuous(breaks = seq(0,1,by = 0.2)) +
  scale_x_continuous(breaks = seq(0,1.4, 0.2)) +
  labs(x = "Frequency of aggression", y = "Presence of new injury",
       color = "Estrous status", linetype = "Estrous status", fill = "Estrous status") +
  theme(text = element_text(family = "Times New Roman"),
        aspect.ratio = 1)  -> p_agginj_all

p_agginj_all

ggsave("figures/p_agginj_all.png", p_agginj_all,
       width = 120, height = 100, units = "mm", dpi = 2000)

実行環境

sessionInfo()
## R version 4.4.0 (2024-04-24 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 22631)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=Japanese_Japan.utf8  LC_CTYPE=Japanese_Japan.utf8   
## [3] LC_MONETARY=Japanese_Japan.utf8 LC_NUMERIC=C                   
## [5] LC_TIME=Japanese_Japan.utf8    
## 
## time zone: Asia/Tokyo
## tzcode source: internal
## 
## attached base packages:
## [1] parallel  stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] systemfonts_1.1.0         extrafont_0.19           
##  [3] ggpattern_1.1.1           ggh4x_0.2.8              
##  [5] ggsignif_0.6.4            lemon_0.4.9              
##  [7] ggsci_3.2.0               kableExtra_1.4.0         
##  [9] patchwork_1.2.0           ggrepel_0.9.6            
## [11] flextable_0.9.6           GGally_2.2.1             
## [13] ggnewscale_0.4.10         viridis_0.6.5            
## [15] viridisLite_0.4.2         bayesplot_1.11.1         
## [17] plotly_4.10.4             ggbeeswarm_0.7.2         
## [19] ggforce_0.4.2             ggtext_0.1.2             
## [21] ggraph_2.2.1              tidygraph_1.3.1          
## [23] igraph_2.0.3              vegan_2.6-6.1            
## [25] lattice_0.22-6            permute_0.9-7            
## [27] aninet_0.2.0.0            asnipe_1.1.17            
## [29] ANTs_0.0.16               sna_2.7-2                
## [31] network_1.18.2            statnet.common_4.9.0     
## [33] DHARMa.helpers_0.0.0.9000 DHARMa_0.4.6             
## [35] ggeffects_1.6.0           cmdstanr_0.8.0           
## [37] rstanarm_2.32.1           rstan_2.32.6             
## [39] StanHeaders_2.32.8        rethinking_2.13.2        
## [41] brms_2.21.0               Rcpp_1.0.12              
## [43] emmeans_1.10.2            scales_1.3.0             
## [45] DT_0.33                   see_0.8.4                
## [47] report_0.5.8              parameters_0.22.0        
## [49] performance_0.12.0        modelbased_0.8.8         
## [51] insight_0.20.1            effectsize_0.8.9         
## [53] datawizard_0.11.0         correlation_0.8.5        
## [55] bayestestR_0.13.2         easystats_0.7.2          
## [57] knitr_1.48                readxl_1.4.3             
## [59] lubridate_1.9.3           forcats_1.0.0            
## [61] stringr_1.5.1             dplyr_1.1.4              
## [63] purrr_1.0.2               readr_2.1.5              
## [65] tidyr_1.3.1               tibble_3.2.1             
## [67] ggplot2_3.5.1             tidyverse_2.0.0          
## 
## loaded via a namespace (and not attached):
##   [1] matrixStats_1.3.0       doParallel_1.0.17       httr_1.4.7             
##   [4] RColorBrewer_1.1-3      numDeriv_2016.8-1.1     tools_4.4.0            
##   [7] backports_1.5.0         sjlabelled_1.2.0        utf8_1.2.4             
##  [10] R6_2.5.1                lazyeval_0.2.2          mgcv_1.9-1             
##  [13] withr_3.0.1             Brobdingnag_1.2-9       gridExtra_2.3          
##  [16] archive_1.1.8           cli_3.6.2               textshaping_0.4.0      
##  [19] shinyjs_2.1.0           officer_0.6.6           sandwich_3.1-0         
##  [22] labeling_0.4.3          sass_0.4.9              mvtnorm_1.2-5          
##  [25] askpass_1.2.0           QuickJSR_1.1.3          gfonts_0.2.0           
##  [28] svglite_2.1.3           rstudioapi_0.16.0       httpcode_0.3.0         
##  [31] generics_0.1.3          shape_1.4.6.1           vroom_1.6.5            
##  [34] gtools_3.9.5            crosstalk_1.2.1         distributional_0.4.0   
##  [37] zip_2.3.1               inline_0.3.19           loo_2.7.0              
##  [40] Matrix_1.7-0            fansi_1.0.6             abind_1.4-8            
##  [43] lifecycle_1.0.4         multcomp_1.4-25         yaml_2.3.8             
##  [46] snakecase_0.11.1        grid_4.4.0              qgam_1.3.4             
##  [49] promises_1.3.0          crayon_1.5.2            miniUI_0.1.1.1         
##  [52] haven_2.5.4             pillar_1.9.0            boot_1.3-30            
##  [55] estimability_1.5.1      shinystan_2.6.0         codetools_0.2-20       
##  [58] glue_1.7.0              V8_4.4.2                fontLiberation_0.1.0   
##  [61] data.table_1.15.4       Rdpack_2.6              vctrs_0.6.5            
##  [64] cellranger_1.1.0        gtable_0.3.5            cachem_1.1.0           
##  [67] xfun_0.44               rbibutils_2.2.16        mime_0.12              
##  [70] coda_0.19-4.1           survival_3.6-4          Kendall_2.2.1          
##  [73] iterators_1.0.14        shinythemes_1.2.0       TH.data_1.1-2          
##  [76] gap_1.5-3               nlme_3.1-164            xts_0.13.2             
##  [79] bit64_4.0.5             fontquiver_0.2.1        threejs_0.3.3          
##  [82] tensorA_0.36.2.1        TMB_1.9.11              bslib_0.7.0            
##  [85] vipor_0.4.7             colorspace_2.1-0        tidyselect_1.2.1       
##  [88] processx_3.8.4          bit_4.0.5               extrafontdb_1.0        
##  [91] compiler_4.4.0          curl_5.2.1              xml2_1.3.6             
##  [94] fontBitstreamVera_0.1.1 colourpicker_1.3.0      posterior_1.5.0        
##  [97] bookdown_0.39           checkmate_2.3.1         dygraphs_1.1.1.6       
## [100] digest_0.6.35           minqa_1.2.7             rmarkdown_2.27         
## [103] htmltools_0.5.8.1       pkgconfig_2.0.3         base64enc_0.1-3        
## [106] lme4_1.1-35.3           highr_0.11              fastmap_1.2.0          
## [109] rlang_1.1.3             htmlwidgets_1.6.4       shiny_1.8.1.1          
## [112] farver_2.1.2            jquerylib_0.1.4         zoo_1.8-12             
## [115] jsonlite_1.8.8          magrittr_2.0.3          munsell_0.5.1          
## [118] gdtools_0.3.7           stringi_1.8.4           MASS_7.3-60.2          
## [121] plyr_1.8.9              pkgbuild_1.4.4          ggstats_0.6.0          
## [124] graphlayouts_1.1.1      splines_4.4.0           gridtext_0.1.5         
## [127] hms_1.1.3               ps_1.7.6                uuid_1.2-0             
## [130] markdown_1.13           reshape2_1.4.4          stats4_4.4.0           
## [133] rstantools_2.4.0        crul_1.4.2              evaluate_1.0.0         
## [136] RcppParallel_5.1.7      foreach_1.5.2           nloptr_2.0.3           
## [139] tzdb_0.4.0              tweenr_2.0.3            httpuv_1.6.15          
## [142] Rttf2pt1_1.3.12         openssl_2.2.0           polyclip_1.10-6        
## [145] gap.datasets_0.0.6      xtable_1.8-4            later_1.3.2            
## [148] ragg_1.3.2              glmmTMB_1.1.9           memoise_2.0.1          
## [151] beeswarm_0.4.0          cluster_2.1.6           timechange_0.3.0       
## [154] bridgesampling_1.1-2

References

Chang, W. (2018). R graphics cookbook: Practical recipes for visualizing data. “O’Reilly Media, Inc.”
Gabry, J., & Mahr, T. (2022). Bayesplot: Plotting for bayesian models. https://mc-stan.org/bayesplot
Hartig, F. (2022). DHARMa: Residual diagnostics for hierarchical (multi-level / mixed) regression models. https://CRAN.R-project.org/package=DHARMa
Rodríguez-Sánchez, F. (2023). DHARMa.helpers: Helper functions to check models not (yet) directly supported by DHARMa.
Wickham, H., & Grolemund, G. (2016). R for data science: Import, tidy, transform, visualize, and model data. “O’Reilly Media, Inc.”
Yamaguchi, T., & Kazahari, N. (2022). Fission–fusion dynamics in a wild group of japanese macaques (macaca fuscata) on kinkazan island caused by the repeated separation of an alpha male being followed by females. Primates, 63(6), 575–582.
松村優哉., 湯谷啓明., 紀ノ定保礼., & 前田和. (2021). RユーザのためのRstudio[実践]入門 tidyverseによるモダンな分析フローの世界 改訂2版. 技術評論社.